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

Optimize SVF calculation, clean-up, bugfixes

File:
1 edited

Legend:

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

    r2906 r2920  
    1616!
    1717! Copyright 2015-2018 Czech Technical University in Prague
     18! Copyright 2015-2018 Institute of Computer Science of the
     19!                     Czech Academy of Sciences, Prague
    1820! Copyright 1997-2018 Leibniz Universitaet Hannover
    1921!------------------------------------------------------------------------------!
     
    2628! -----------------
    2729! $Id$
     30! Remove unused pcbl, npcbl from ONLY list
     31! moh.hefny:
     32! Fixed bugs introduced by new structures and by moving radiation interaction
     33! into radiation_model_mod.f90.
     34! Bugfix: usm data output 3D didn't respect directions
     35!
     36! 2906 2018-03-19 08:56:40Z Giersch
    2837! Local variable ids has to be initialized with a value of -1 in
    2938! usm_average_3d_data
     
    217226!> Further work:
    218227!> -------------
    219 !> 1. Reduce number of shape view factors by merging factors for distant surfaces
    220 !>    under shallow angles. Idea: Iteratively select the smallest shape view
    221 !>    factor by value (among all sources and targets) which has a similarly
    222 !>    oriented source neighbor (or near enough) SVF and merge them by adding
    223 !>    value of the smaller SVF to the larger one and deleting the smaller one.
    224 !>    This will allow for better scaling at higher resolutions.
    225 !>
    226 !> 2. Remove global arrays surfouts, surfoutl and only keep track of radiosity
     228!> 1. Remove global arrays surfouts, surfoutl and only keep track of radiosity
    227229!>    from surfaces that are visible from local surfaces (i.e. there is a SVF
    228230!>    where target is local). To do that, radiosity will be exchanged after each
    229231!>    reflection step using MPI_Alltoall instead of current MPI_Allgather.
    230232!>
    231 !> 3. Temporarily large values of surface heat flux can be observed, up to
     233!> 2. Temporarily large values of surface heat flux can be observed, up to
    232234!>    1.2 Km/s, which seem to be not realistic.
    233235!>
     
    248250#if ! defined( __nopointer )
    249251    USE arrays_3d,                                                             &
    250         ONLY:  zu, pt, pt_1, pt_2, p, u, v, w, hyp, tend
     252        ONLY:  hyp, zu, pt, pt_1, pt_2, p, u, v, w, hyp, tend
    251253#endif
    252254
     
    272274
    273275    USE date_and_time_mod,                                                     &
    274         ONLY:  d_seconds_year, day_of_year_init, time_utc_init
     276        ONLY:  time_utc_init
    275277
    276278    USE grid_variables,                                                        &
     
    288290   
    289291    USE plant_canopy_model_mod,                                                &
    290         ONLY:  pc_heating_rate, usm_lad_rma
     292        ONLY:  pc_heating_rate
    291293   
    292294    USE radiation_model_mod,                                                   &
    293         ONLY:  albedo_type, radiation, calc_zenith, zenith,                    &
     295        ONLY:  albedo_type, radiation_interaction, calc_zenith, zenith,                    &
    294296               rad_sw_in, rad_lw_in, rad_sw_out, rad_lw_out,                   &
    295297               sigma_sb, solar_constant, sun_direction, sun_dir_lat,           &
     
    297299               force_radiation_call, surfinsw, surfinlw, surfinswdir,          &
    298300               surfinswdif, surfoutsw, surfoutlw, surfins,nsvfl, svf, svfsurf, &
    299                surfinl, surfinlwdif, energy_balance_surf_h,                    &
    300                energy_balance_surf_v, rad_sw_in_dir, rad_sw_in_diff,           &
     301               surfinl, surfinlwdif, rad_sw_in_dir, rad_sw_in_diff,            &
    301302               rad_lw_in_diff, surfouts, surfoutl, surfoutsl, surfoutll, surf, &
    302                surfl, nsurfl, nsurfs, surfstart, pcbinsw, pcbinlw, pcbl, npcbl, startenergy,      &
    303                endenergy, nenergy, iup_u, inorth_u, isouth_u, ieast_u, iwest_u, iup_l,  &
    304                inorth_l, isouth_l, ieast_l, iwest_l, startsky, endsky,id,      &
    305                iz, iy, ix, idir, jdir, kdir, startborder, endborder, nsurf_type, nzub, nzut,   &
    306                isky, inorth_b,idown_a, isouth_b, ieast_b, iwest_b, nzu, pch, nsurf,  &
    307                iup_a, inorth_a, isouth_a, ieast_a, iwest_a, idsvf, ndsvf,      &
    308                idcsf, ndcsf, kdcsf, pct, startland, endland, startwall, endwall
     303               surfl, nsurfl, nsurfs, surfstart, pcbinsw, pcbinlw,             &
     304               iup_u, inorth_u, isouth_u, ieast_u, iwest_u, iup_l,             &
     305               inorth_l, isouth_l, ieast_l, iwest_l, id,                       &
     306               iz, iy, ix, idir, jdir, kdir,  nsurf_type, nzub, nzut,          &
     307               nzu, pch, nsurf, idsvf, ndsvf,                                  &
     308               iup_a, idown_a, inorth_a, isouth_a, ieast_a, iwest_a,           &
     309               idcsf, ndcsf, kdcsf, pct,                                       &
     310               startland, endland, startwall, endwall, skyvf, skyvft
    309311
    310312    USE statistics,                                                            &
     
    322324    LOGICAL                                        ::  force_radiation_call_l = .FALSE.   !< flag parameter for unscheduled radiation model calls
    323325    LOGICAL                                        ::  indoor_model = .FALSE.              !< whether to use the indoor model
    324 
    325    
     326    LOGICAL                                        ::  read_wall_temp_3d = .FALSE.
     327
     328
    326329    INTEGER(iwp)                                   ::  building_type = 1                  !< default building type (preleminary setting)
    327330    INTEGER(iwp)                                   ::  land_category = 2                  !< default category for land surface
    328331    INTEGER(iwp)                                   ::  wall_category = 2                  !< default category for wall surface over pedestrian zone
    329     INTEGER(iwp)                                   ::  pedestrant_category = 2            !< default category for wall surface in pedestrian zone
     332    INTEGER(iwp)                                   ::  pedestrian_category = 2            !< default category for wall surface in pedestrian zone
    330333    INTEGER(iwp)                                   ::  roof_category = 2                  !< default category for root surface
    331 
     334    REAL(wp)                                       ::  roughness_concrete = 0.001_wp      !< roughness length of average concrete surface
    332335!
    333336!-- Indices of input attributes for (above) ground floor level
     
    380383    REAL(wp)  ::  roof_height_limit = 4._wp          !< height for distinguish between land surfaces and roofs
    381384    REAL(wp)  ::  ground_floor_level = 4.0_wp        !< default ground floor level
    382     REAL(wp)  ::  ra_horiz_coef = 5.0_wp             !< mysterious coefficient for correction of overestimation
    383                                                                                           !< of r_a for horizontal surfaces -> TODO
     385
    384386
    385387    CHARACTER(37), DIMENSION(0:6), PARAMETER :: building_type_name = (/     &
     
    524526!-- anthropogenic heat sources
    525527!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    526     REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  aheat             !< daily average of anthropogenic heat (W/m2)
    527     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  aheatprof         !< diurnal profile of anthropogenic heat
     528    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE        ::  aheat             !< daily average of anthropogenic heat (W/m2)
     529    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  aheatprof         !< diurnal profiles of anthropogenic heat for particular layers
     530    INTEGER(wp)                                    ::  naheatlayers = 1  !< number of layers of anthropogenic heat
    528531
    529532!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     
    547550!    REAL(wp), DIMENSION(nzb_wall:nzt_wall)         :: zwn_default_green = (/0.33_wp, 0.66_wp, 1.0_wp /)
    548551
    549                                                                        
    550     REAL(wp)                                       ::   wall_inner_temperature = 296.0_wp    !< temperature of the inner wall surface (~23 degrees C) (K)
    551     REAL(wp)                                       ::   roof_inner_temperature = 296.0_wp    !< temperature of the inner roof surface (~23 degrees C) (K)
    552     REAL(wp)                                       ::   soil_inner_temperature = 283.0_wp    !< temperature of the deep soil (~10 degrees C) (K)
    553     REAL(wp)                                       ::   window_inner_temperature = 296.0_wp  !< temperature of the inner window surface (~23 degrees C) (K)
     552
     553    REAL(wp)                                       :: wall_inner_temperature = 295.0_wp    !< temperature of the inner wall surface (~22 degrees C) (K)
     554    REAL(wp)                                       :: roof_inner_temperature = 295.0_wp    !< temperature of the inner roof surface (~22 degrees C) (K)
     555    REAL(wp)                                       :: soil_inner_temperature = 288.0_wp    !< temperature of the deep soil (~15 degrees C) (K)
     556    REAL(wp)                                       :: window_inner_temperature = 295.0_wp  !< temperature of the inner window surface (~22 degrees C) (K)
    554557
    555558!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     
    674677!-- albedo, emissivity, lambda_surf, roughness, thickness, volumetric heat capacity, thermal conductivity
    675678    INTEGER(iwp)                                   :: n_surface_types      !< number of the wall type categories
    676     INTEGER(iwp), PARAMETER                        :: n_surface_params = 8 !< number of parameters for each type of the wall
     679    INTEGER(iwp), PARAMETER                        :: n_surface_params = 9 !< number of parameters for each type of the wall
    677680    INTEGER(iwp), PARAMETER                        :: ialbedo  = 1         !< albedo of the surface
    678681    INTEGER(iwp), PARAMETER                        :: iemiss   = 2         !< emissivity of the surface
    679     INTEGER(iwp), PARAMETER                        :: ilambdas = 3         !< heat conductivity λS between air and surface ( W m−2 K−1 )
    680     INTEGER(iwp), PARAMETER                        :: irough   = 4         !< roughness relative to concrete
    681     INTEGER(iwp), PARAMETER                        :: icsurf   = 5         !< Surface skin layer heat capacity (J m−2 K−1 )
    682     INTEGER(iwp), PARAMETER                        :: ithick   = 6         !< thickness of the surface (wall, roof, land)  ( m )
    683     INTEGER(iwp), PARAMETER                        :: irhoC    = 7         !< volumetric heat capacity rho*C of the material ( J m−3 K−1 )
    684     INTEGER(iwp), PARAMETER                        :: ilambdah = 8         !< thermal conductivity λH of the wall (W m−1 K−1 )
     682    INTEGER(iwp), PARAMETER                        :: ilambdas = 3         !< heat conductivity lambda S between surface and material ( W m-2 K-1 )
     683    INTEGER(iwp), PARAMETER                        :: irough   = 4         !< roughness length z0 for movements
     684    INTEGER(iwp), PARAMETER                        :: iroughh  = 5         !< roughness length z0h for scalars (heat, humidity,...)
     685    INTEGER(iwp), PARAMETER                        :: icsurf   = 6         !< Surface skin layer heat capacity (J m-2 K-1 )
     686    INTEGER(iwp), PARAMETER                        :: ithick   = 7         !< thickness of the surface (wall, roof, land)  ( m )
     687    INTEGER(iwp), PARAMETER                        :: irhoC    = 8         !< volumetric heat capacity rho*C of the material ( J m-3 K-1 )
     688    INTEGER(iwp), PARAMETER                        :: ilambdah = 9         !< thermal conductivity lambda H of the wall (W m-1 K-1 )
    685689    CHARACTER(12), DIMENSION(:), ALLOCATABLE       :: surface_type_names   !< names of wall types (used only for reports)
    686690    INTEGER(iwp), DIMENSION(:), ALLOCATABLE        :: surface_type_codes   !< codes of wall types
     
    760764!-- Public functions
    761765    PUBLIC usm_boundary_condition, usm_check_parameters, usm_init_urban_surface,&
    762            usm_rrd_local,                                                       &
    763            usm_surface_energy_balance, usm_material_heat_model,                 &
    764            usm_swap_timelevel, usm_check_data_output, usm_average_3d_data,      &
    765            usm_data_output_3d, usm_define_netcdf_grid, usm_parin,               &
     766           usm_rrd_local,                                                      &
     767           usm_surface_energy_balance, usm_material_heat_model,                &
     768           usm_swap_timelevel, usm_check_data_output, usm_average_3d_data,     &
     769           usm_data_output_3d, usm_define_netcdf_grid, usm_parin,              &
    766770           usm_wrd_local, usm_allocate_surface
    767771
    768772!-- Public parameters, constants and initial values
    769     PUBLIC usm_anthropogenic_heat, usm_material_model, ra_horiz_coef,           &
     773    PUBLIC usm_anthropogenic_heat, usm_material_model,                          &
    770774           usm_green_heat_model, usm_temperature_near_surface
    771775
     
    11751179        CHARACTER (len=*), INTENT(IN) :: variable
    11761180 
    1177         INTEGER(iwp)                                       :: i, j, k, l, m, ids, iwl,istat
     1181        INTEGER(iwp)                                       :: i, j, k, l, m, ids, idsint, iwl, istat
    11781182        CHARACTER (len=varnamelength)                      :: var, surfid
    11791183        INTEGER(iwp), PARAMETER                            :: nd = 5
    11801184        CHARACTER(len=6), DIMENSION(0:nd-1), PARAMETER     :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /)
     1185        INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER         :: dirint = (/ iup_u, isouth_u, inorth_u, iwest_u, ieast_u /)
    11811186
    11821187!--     find the real name of the variable
     
    11881193            IF ( var(k-j+1:k) == dirname(i) )  THEN
    11891194                ids = i
     1195                idsint = dirint(ids)
    11901196                var = var(:k-j)
    11911197                EXIT
     
    12451251                CASE ( 'usm_rad_insw' )
    12461252!--                 array of sw radiation falling to surface after i-th reflection
    1247                     IF ( .NOT.  ALLOCATED(surf_usm_h%surfinsw_av) )  THEN
    1248                         ALLOCATE( surf_usm_h%surfinsw_av(1:surf_usm_h%ns) )
    1249                         surf_usm_h%surfinsw_av = 0.0_wp
     1253                    IF ( .NOT.  ALLOCATED(surfinsw_av) )  THEN
     1254                        ALLOCATE( surfinsw_av(nsurfl) )
     1255                        surfinsw_av = 0.0_wp
    12501256                    ENDIF
    1251                     DO  l = 0, 3
    1252                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%surfinsw_av) )  THEN
    1253                            ALLOCATE( surf_usm_v(l)%surfinsw_av(1:surf_usm_v(l)%ns) )
    1254                            surf_usm_v(l)%surfinsw_av = 0.0_wp
    1255                        ENDIF
    1256                     ENDDO
    1257                                    
     1257
    12581258                CASE ( 'usm_rad_inlw' )
    12591259!--                 array of lw radiation falling to surface after i-th reflection
    1260                     IF ( .NOT.  ALLOCATED(surf_usm_h%surfinlw_av) )  THEN
    1261                         ALLOCATE( surf_usm_h%surfinlw_av(1:surf_usm_h%ns) )
    1262                         surf_usm_h%surfinlw_av = 0.0_wp
     1260                    IF ( .NOT.  ALLOCATED(surfinlw_av) )  THEN
     1261                        ALLOCATE( surfinlw_av(nsurfl) )
     1262                        surfinlw_av = 0.0_wp
    12631263                    ENDIF
    1264                     DO  l = 0, 3
    1265                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%surfinlw_av) )  THEN
    1266                            ALLOCATE( surf_usm_v(l)%surfinlw_av(1:surf_usm_v(l)%ns) )
    1267                            surf_usm_v(l)%surfinlw_av = 0.0_wp
    1268                        ENDIF
    1269                     ENDDO
    12701264
    12711265                CASE ( 'usm_rad_inswdir' )
    12721266!--                 array of direct sw radiation falling to surface from sun
    12731267                    IF ( .NOT.  ALLOCATED(surfinswdir_av) )  THEN
    1274                         ALLOCATE( surfinswdir_av(startenergy:endenergy) )
     1268                        ALLOCATE( surfinswdir_av(nsurfl) )
    12751269                        surfinswdir_av = 0.0_wp
    12761270                    ENDIF
     
    12791273!--                 array of difusion sw radiation falling to surface from sky and borders of the domain
    12801274                    IF ( .NOT.  ALLOCATED(surfinswdif_av) )  THEN
    1281                         ALLOCATE( surfinswdif_av(startenergy:endenergy) )
     1275                        ALLOCATE( surfinswdif_av(nsurfl) )
    12821276                        surfinswdif_av = 0.0_wp
    12831277                    ENDIF
     
    12861280!--                 array of sw radiation falling to surface from reflections
    12871281                    IF ( .NOT.  ALLOCATED(surfinswref_av) )  THEN
    1288                         ALLOCATE( surfinswref_av(startenergy:endenergy) )
     1282                        ALLOCATE( surfinswref_av(nsurfl) )
    12891283                        surfinswref_av = 0.0_wp
    12901284                    ENDIF
     
    12931287!--                 array of sw radiation falling to surface after i-th reflection
    12941288                   IF ( .NOT.  ALLOCATED(surfinlwdif_av) )  THEN
    1295                         ALLOCATE( surfinlwdif_av(startenergy:endenergy) )
     1289                        ALLOCATE( surfinlwdif_av(nsurfl) )
    12961290                        surfinlwdif_av = 0.0_wp
    12971291                    ENDIF
     
    13001294!--                 array of lw radiation falling to surface from reflections
    13011295                    IF ( .NOT.  ALLOCATED(surfinlwref_av) )  THEN
    1302                         ALLOCATE( surfinlwref_av(startenergy:endenergy) )
     1296                        ALLOCATE( surfinlwref_av(nsurfl) )
    13031297                        surfinlwref_av = 0.0_wp
    13041298                    ENDIF
     
    13071301!--                 array of sw radiation emitted from surface after i-th reflection
    13081302                    IF ( .NOT.  ALLOCATED(surfoutsw_av) )  THEN
    1309                         ALLOCATE( surfoutsw_av(startenergy:endenergy) )
     1303                        ALLOCATE( surfoutsw_av(nsurfl) )
    13101304                        surfoutsw_av = 0.0_wp
    13111305                    ENDIF
     
    13141308!--                 array of lw radiation emitted from surface after i-th reflection
    13151309                    IF ( .NOT.  ALLOCATED(surfoutlw_av) )  THEN
    1316                         ALLOCATE( surfoutlw_av(startenergy:endenergy) )
     1310                        ALLOCATE( surfoutlw_av(nsurfl) )
    13171311                        surfoutlw_av = 0.0_wp
    13181312                    ENDIF
     
    13201314!--                 array of residua of sw radiation absorbed in surface after last reflection
    13211315                    IF ( .NOT.  ALLOCATED(surfins_av) )  THEN
    1322                         ALLOCATE( surfins_av(startenergy:endenergy) )
     1316                        ALLOCATE( surfins_av(nsurfl) )
    13231317                        surfins_av = 0.0_wp
    13241318                    ENDIF
     
    13271321!--                 array of residua of lw radiation absorbed in surface after last reflection
    13281322                    IF ( .NOT.  ALLOCATED(surfinl_av) )  THEN
    1329                         ALLOCATE( surfinl_av(startenergy:endenergy) )
     1323                        ALLOCATE( surfinl_av(nsurfl) )
    13301324                        surfinl_av = 0.0_wp
    13311325                    ENDIF
     
    15441538                CASE ( 'usm_rad_insw' )
    15451539!--                 array of sw radiation falling to surface after i-th reflection
    1546                     DO l = startenergy, endenergy
    1547                         IF ( surfl(id,l) == ids )  THEN
     1540                    DO l = 1, nsurfl
     1541                        IF ( surfl(id,l) == idsint )  THEN
    15481542                            surfinsw_av(l) = surfinsw_av(l) + surfinsw(l)
    15491543                        ENDIF
     
    15521546                CASE ( 'usm_rad_inlw' )
    15531547!--                 array of lw radiation falling to surface after i-th reflection
    1554                     DO l = startenergy, endenergy
    1555                         IF ( surfl(id,l) == ids )  THEN
     1548                    DO l = 1, nsurfl
     1549                        IF ( surfl(id,l) == idsint )  THEN
    15561550                            surfinlw_av(l) = surfinlw_av(l) + surfinlw(l)
    15571551                        ENDIF
     
    15601554                CASE ( 'usm_rad_inswdir' )
    15611555!--                 array of direct sw radiation falling to surface from sun
    1562                     DO l = startenergy, endenergy
    1563                         IF ( surfl(id,l) == ids )  THEN
     1556                    DO l = 1, nsurfl
     1557                        IF ( surfl(id,l) == idsint )  THEN
    15641558                            surfinswdir_av(l) = surfinswdir_av(l) + surfinswdir(l)
    15651559                        ENDIF
     
    15681562                CASE ( 'usm_rad_inswdif' )
    15691563!--                 array of difusion sw radiation falling to surface from sky and borders of the domain
    1570                     DO l = startenergy, endenergy
    1571                         IF ( surfl(id,l) == ids )  THEN
     1564                    DO l = 1, nsurfl
     1565                        IF ( surfl(id,l) == idsint )  THEN
    15721566                            surfinswdif_av(l) = surfinswdif_av(l) + surfinswdif(l)
    15731567                        ENDIF
     
    15761570                CASE ( 'usm_rad_inswref' )
    15771571!--                 array of sw radiation falling to surface from reflections
    1578                     DO l = startenergy, endenergy
    1579                         IF ( surfl(id,l) == ids )  THEN
     1572                    DO l = 1, nsurfl
     1573                        IF ( surfl(id,l) == idsint )  THEN
    15801574                            surfinswref_av(l) = surfinswref_av(l) + surfinsw(l) - &
    15811575                                                surfinswdir(l) - surfinswdif(l)
     
    15861580                CASE ( 'usm_rad_inlwdif' )
    15871581!--                 array of sw radiation falling to surface after i-th reflection
    1588                     DO l = startenergy, endenergy
    1589                         IF ( surfl(id,l) == ids )  THEN
    1590                             surfinswref_av(l) = surfinswref_av(l) + surfinsw(l) - &
    1591                                                 surfinswdir(l) - surfinswdif(l)
     1582                    DO l = 1, nsurfl
     1583                        IF ( surfl(id,l) == idsint )  THEN
     1584                            surfinlwdif_av(l) = surfinlwdif_av(l) + surfinlwdif(l)
    15921585                        ENDIF
    15931586                    ENDDO
     
    15951588                CASE ( 'usm_rad_inlwref' )
    15961589!--                 array of lw radiation falling to surface from reflections
    1597                     DO l = startenergy, endenergy
    1598                         IF ( surfl(id,l) == ids )  THEN
    1599                             surfinlwdif_av(l) = surfinlwdif_av(l) + surfinlwdif(l)
     1590                    DO l = 1, nsurfl
     1591                        IF ( surfl(id,l) == idsint )  THEN
     1592                            surfinlwref_av(l) = surfinlwref_av(l) + &
     1593                                                surfinlw(l) - surfinlwdif(l)
    16001594                        ENDIF
    16011595                    ENDDO
     
    16031597                CASE ( 'usm_rad_outsw' )
    16041598!--                 array of sw radiation emitted from surface after i-th reflection
    1605                     DO l = startenergy, endenergy
    1606                         IF ( surfl(id,l) == ids )  THEN
    1607                             surfinlwref_av(l) = surfinlwref_av(l) + &
    1608                                                 surfinlw(l) - surfinlwdif(l)
     1599                    DO l = 1, nsurfl
     1600                        IF ( surfl(id,l) == idsint )  THEN
     1601                            surfoutsw_av(l) = surfoutsw_av(l) + surfoutsw(l)
    16091602                        ENDIF
    16101603                    ENDDO
     
    16121605                CASE ( 'usm_rad_outlw' )
    16131606!--                 array of lw radiation emitted from surface after i-th reflection
    1614                     DO l = startenergy, endenergy
    1615                         IF ( surfl(id,l) == ids )  THEN
    1616                             surfoutsw_av(l) = surfoutsw_av(l) + surfoutsw(l)
     1607                    DO l = 1, nsurfl
     1608                        IF ( surfl(id,l) == idsint )  THEN
     1609                            surfoutlw_av(l) = surfoutlw_av(l) + surfoutlw(l)
    16171610                        ENDIF
    16181611                    ENDDO
     
    16201613                CASE ( 'usm_rad_ressw' )
    16211614!--                 array of residua of sw radiation absorbed in surface after last reflection
    1622                     DO l = startenergy, endenergy
    1623                         IF ( surfl(id,l) == ids )  THEN
    1624                             surfoutlw_av(l) = surfoutlw_av(l) + surfoutlw(l)
     1615                    DO l = 1, nsurfl
     1616                        IF ( surfl(id,l) == idsint )  THEN
     1617                            surfins_av(l) = surfins_av(l) + surfins(l)
    16251618                        ENDIF
    16261619                    ENDDO
     
    16281621                CASE ( 'usm_rad_reslw' )
    16291622!--                 array of residua of lw radiation absorbed in surface after last reflection
    1630                     DO l = startenergy, endenergy
    1631                         IF ( surfl(id,l) == ids )  THEN
    1632                             surfins_av(l) = surfins_av(l) + surfins(l)
     1623                    DO l = 1, nsurfl
     1624                        IF ( surfl(id,l) == idsint )  THEN
     1625                            surfinl_av(l) = surfinl_av(l) + surfinl(l)
    16331626                        ENDIF
    16341627                    ENDDO
     
    18711864                CASE ( 'usm_rad_insw' )
    18721865!--                 array of sw radiation falling to surface after i-th reflection
    1873                     DO l = startenergy, endenergy
    1874                         IF ( surfl(id,l) == ids )  THEN
     1866                    DO l = 1, nsurfl
     1867                        IF ( surfl(id,l) == idsint )  THEN
    18751868                            surfinsw_av(l) = surfinsw_av(l) / REAL( average_count_3d, kind=wp )
    18761869                        ENDIF
     
    18791872                CASE ( 'usm_rad_inlw' )
    18801873!--                 array of lw radiation falling to surface after i-th reflection
    1881                     DO l = startenergy, endenergy
    1882                         IF ( surfl(id,l) == ids )  THEN
     1874                    DO l = 1, nsurfl
     1875                        IF ( surfl(id,l) == idsint )  THEN
    18831876                            surfinlw_av(l) = surfinlw_av(l) / REAL( average_count_3d, kind=wp )
    18841877                        ENDIF
     
    18871880                CASE ( 'usm_rad_inswdir' )
    18881881!--                 array of direct sw radiation falling to surface from sun
    1889                     DO l = startenergy, endenergy
    1890                         IF ( surfl(id,l) == ids )  THEN
     1882                    DO l = 1, nsurfl
     1883                        IF ( surfl(id,l) == idsint )  THEN
    18911884                            surfinswdir_av(l) = surfinswdir_av(l) / REAL( average_count_3d, kind=wp )
    18921885                        ENDIF
     
    18951888                CASE ( 'usm_rad_inswdif' )
    18961889!--                 array of difusion sw radiation falling to surface from sky and borders of the domain
    1897                     DO l = startenergy, endenergy
    1898                         IF ( surfl(id,l) == ids )  THEN
     1890                    DO l = 1, nsurfl
     1891                        IF ( surfl(id,l) == idsint )  THEN
    18991892                            surfinswdif_av(l) = surfinswdif_av(l) / REAL( average_count_3d, kind=wp )
    19001893                        ENDIF
     
    19031896                CASE ( 'usm_rad_inswref' )
    19041897!--                 array of sw radiation falling to surface from reflections
    1905                     DO l = startenergy, endenergy
    1906                         IF ( surfl(id,l) == ids )  THEN
     1898                    DO l = 1, nsurfl
     1899                        IF ( surfl(id,l) == idsint )  THEN
    19071900                            surfinswref_av(l) = surfinswref_av(l) / REAL( average_count_3d, kind=wp )
    19081901                        ENDIF
     
    19111904                CASE ( 'usm_rad_inlwdif' )
    19121905!--                 array of sw radiation falling to surface after i-th reflection
    1913                     DO l = startenergy, endenergy
    1914                         IF ( surfl(id,l) == ids )  THEN
     1906                    DO l = 1, nsurfl
     1907                        IF ( surfl(id,l) == idsint )  THEN
    19151908                            surfinlwdif_av(l) = surfinlwdif_av(l) / REAL( average_count_3d, kind=wp )
    19161909                        ENDIF
     
    19191912                CASE ( 'usm_rad_inlwref' )
    19201913!--                 array of lw radiation falling to surface from reflections
    1921                     DO l = startenergy, endenergy
    1922                         IF ( surfl(id,l) == ids )  THEN
     1914                    DO l = 1, nsurfl
     1915                        IF ( surfl(id,l) == idsint )  THEN
    19231916                            surfinlwref_av(l) = surfinlwref_av(l) / REAL( average_count_3d, kind=wp )
    19241917                        ENDIF
     
    19271920                CASE ( 'usm_rad_outsw' )
    19281921!--                 array of sw radiation emitted from surface after i-th reflection
    1929                     DO l = startenergy, endenergy
    1930                         IF ( surfl(id,l) == ids )  THEN
     1922                    DO l = 1, nsurfl
     1923                        IF ( surfl(id,l) == idsint )  THEN
    19311924                            surfoutsw_av(l) = surfoutsw_av(l) / REAL( average_count_3d, kind=wp )
    19321925                        ENDIF
     
    19351928                CASE ( 'usm_rad_outlw' )
    19361929!--                 array of lw radiation emitted from surface after i-th reflection
    1937                     DO l = startenergy, endenergy
    1938                         IF ( surfl(id,l) == ids )  THEN
     1930                    DO l = 1, nsurfl
     1931                        IF ( surfl(id,l) == idsint )  THEN
    19391932                            surfoutlw_av(l) = surfoutlw_av(l) / REAL( average_count_3d, kind=wp )
    19401933                        ENDIF
     
    19431936                CASE ( 'usm_rad_ressw' )
    19441937!--                 array of residua of sw radiation absorbed in surface after last reflection
    1945                     DO l = startenergy, endenergy
    1946                         IF ( surfl(id,l) == ids )  THEN
     1938                    DO l = 1, nsurfl
     1939                        IF ( surfl(id,l) == idsint )  THEN
    19471940                            surfins_av(l) = surfins_av(l) / REAL( average_count_3d, kind=wp )
    19481941                        ENDIF
     
    19511944                CASE ( 'usm_rad_reslw' )
    19521945!--                 array of residua of lw radiation absorbed in surface after last reflection
    1953                     DO l = startenergy, endenergy
    1954                         IF ( surfl(id,l) == ids )  THEN
     1946                    DO l = 1, nsurfl
     1947                        IF ( surfl(id,l) == idsint )  THEN
    19551948                            surfinl_av(l) = surfinl_av(l) / REAL( average_count_3d, kind=wp )
    19561949                        ENDIF
     
    22452238             var(1:10) == 'usm_iwghf_' .OR. var(1:17) == 'usm_iwghf_window_' )  THEN
    22462239            unit = 'W/m2'
    2247         ELSE IF ( var(1:10) == 'usm_t_surf'  .OR.  var(1:10) == 'usm_t_wall' .OR.         &
     2240        ELSE IF ( var(1:10) == 'usm_t_surf'   .OR.  var(1:10) == 'usm_t_wall' .OR.         &
    22482241                  var(1:12) == 'usm_t_window' .OR. var(1:17) == 'usm_t_surf_window' .OR.  &
    22492242                  var(1:16) == 'usm_t_surf_green'  .OR.                                   &
     
    22532246        ELSE IF ( var(1:9) == 'usm_surfz'  .OR.  var(1:7) == 'usm_svf'  .OR.              &
    22542247                  var(1:7) == 'usm_dif'  .OR.  var(1:11) == 'usm_surfcat'  .OR.           &
    2255                   var(1:11) == 'usm_surfalb'  .OR.  var(1:12) == 'usm_surfemis')  THEN
     2248                  var(1:11) == 'usm_surfalb'  .OR.  var(1:12) == 'usm_surfemis'  .OR.     &
     2249                  var(1:9) == 'usm_skyvf' .OR. var(1:9) == 'usm_skyvft' )  THEN
    22562250            unit = '1'
    22572251        ELSE
     
    23022296          CALL message( 'check_parameters', 'PA0592', 1, 2, 0, 6, 0 )
    23032297       ENDIF
     2298!
     2299!--    naheatlayers
     2300       IF ( naheatlayers > nzt )  THEN
     2301          message_string = 'number of anthropogenic heat layers '//            &
     2302                           '"naheatlayers" can not be larger than'//           &
     2303                           ' number of domain layers "nzt"'
     2304          CALL message( 'check_parameters', 'PA0593', 1, 2, 0, 6, 0 )
     2305       ENDIF
    23042306
    23052307    END SUBROUTINE usm_check_parameters
     
    23312333        INTEGER(iwp), PARAMETER                                :: nd = 5
    23322334        CHARACTER(len=6), DIMENSION(0:nd-1), PARAMETER         :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /)
    2333         INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER             :: dirint = (/ iup_u, isouth_u, inorth_u, iwest_u, ieast_u /)
     2335        INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER             :: dirint =  (/    iup_u, isouth_u, inorth_u,  iwest_u,  ieast_u /)
     2336        INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER             :: diridx =  (/       -1,        1,        0,        3,        2 /)
     2337                                                                     !< index for surf_*_v: 0:3 = (North, South, East, West)
    23342338        INTEGER(iwp), DIMENSION(0:nd-1)                        :: dirstart
    23352339        INTEGER(iwp), DIMENSION(0:nd-1)                        :: dirend
    2336         INTEGER(iwp)                                           :: ids,isurf,isvf,isurfs,isurflt
     2340        INTEGER(iwp)                                           :: ids,idsint,idsidx,isurf,isvf,isurfs,isurflt
    23372341        INTEGER(iwp)                                           :: is,js,ks,i,j,k,iwl,istat, l, m
    23382342        INTEGER(iwp)                                           ::  k_topo    !< topography top index
     
    23512355            IF ( var(k-j+1:k) == dirname(i) )  THEN
    23522356                ids = i
     2357                idsint = dirint(ids)
     2358                idsidx = diridx(ids)
    23532359                var = var(:k-j)
    23542360                EXIT
     
    24002406          CASE ( 'usm_surfz' )
    24012407!--           array of lw radiation falling to local surface after i-th reflection
    2402               DO  m = 1, surf_usm_h%ns
    2403                  i = surf_usm_h%i(m)
    2404                  j = surf_usm_h%j(m)
    2405                  k = surf_usm_h%k(m)
    2406                  temp_pf(0,j,i) = MAX( temp_pf(0,j,i), REAL( k, kind=wp) )
    2407               ENDDO
    2408               DO  l = 0, 3
     2408              IF ( idsint == iup_u )  THEN
     2409                 DO  m = 1, surf_usm_h%ns
     2410                    i = surf_usm_h%i(m)
     2411                    j = surf_usm_h%j(m)
     2412                    k = surf_usm_h%k(m)
     2413                    temp_pf(0,j,i) = MAX( temp_pf(0,j,i), REAL( k, kind=wp) )
     2414                 ENDDO
     2415              ELSE
     2416                 l = idsidx
    24092417                 DO  m = 1, surf_usm_v(l)%ns
    24102418                    i = surf_usm_v(l)%i(m)
     
    24132421                    temp_pf(0,j,i) = MAX( temp_pf(0,j,i), REAL( k, kind=wp) + 1.0_wp )
    24142422                 ENDDO
    2415               ENDDO
     2423              ENDIF
    24162424
    24172425          CASE ( 'usm_surfcat' )
    24182426!--           surface category
    2419               DO  m = 1, surf_usm_h%ns
    2420                  i = surf_usm_h%i(m)
    2421                  j = surf_usm_h%j(m)
    2422                  k = surf_usm_h%k(m)
    2423                  temp_pf(k,j,i) = surf_usm_h%surface_types(m)
    2424               ENDDO
    2425               DO  l = 0, 3
     2427              IF ( idsint == iup_u )  THEN
     2428                 DO  m = 1, surf_usm_h%ns
     2429                    i = surf_usm_h%i(m)
     2430                    j = surf_usm_h%j(m)
     2431                    k = surf_usm_h%k(m)
     2432                    temp_pf(k,j,i) = surf_usm_h%surface_types(m)
     2433                 ENDDO
     2434              ELSE
     2435                 l = idsidx
    24262436                 DO  m = 1, surf_usm_v(l)%ns
    24272437                    i = surf_usm_v(l)%i(m)
     
    24302440                    temp_pf(k,j,i) = surf_usm_v(l)%surface_types(m)
    24312441                 ENDDO
    2432               ENDDO
     2442              ENDIF
    24332443             
    24342444          CASE ( 'usm_surfalb' )
    24352445!--           surface albedo, weighted average
    2436               DO  m = 1, surf_usm_h%ns
    2437                  i = surf_usm_h%i(m)
    2438                  j = surf_usm_h%j(m)
    2439                  k = surf_usm_h%k(m)
    2440                  temp_pf(k,j,i) = surf_usm_h%frac(0,m)     *                   &
    2441                                   surf_usm_h%albedo(0,m) +                     &
    2442                                   surf_usm_h%frac(1,m)     *                   &
    2443                                   surf_usm_h%albedo(1,m) +                     &
    2444                                   surf_usm_h%frac(2,m)     *                   &
    2445                                   surf_usm_h%albedo(2,m)
    2446               ENDDO
    2447               DO  l = 0, 3
     2446              IF ( idsint == iup_u )  THEN
     2447                 DO  m = 1, surf_usm_h%ns
     2448                    i = surf_usm_h%i(m)
     2449                    j = surf_usm_h%j(m)
     2450                    k = surf_usm_h%k(m)
     2451                    temp_pf(k,j,i) = surf_usm_h%frac(0,m)     *                   &
     2452                                     surf_usm_h%albedo(0,m) +                     &
     2453                                     surf_usm_h%frac(1,m)     *                   &
     2454                                     surf_usm_h%albedo(1,m) +                     &
     2455                                     surf_usm_h%frac(2,m)     *                   &
     2456                                     surf_usm_h%albedo(2,m)
     2457                 ENDDO
     2458              ELSE
     2459                 l = idsidx
    24482460                 DO  m = 1, surf_usm_v(l)%ns
    24492461                    i = surf_usm_v(l)%i(m)
     
    24572469                                     surf_usm_v(l)%albedo(2,m)
    24582470                 ENDDO
    2459               ENDDO
     2471              ENDIF
    24602472             
    24612473          CASE ( 'usm_surfemis' )
    24622474!--           surface emissivity, weighted average
    2463               DO  m = 1, surf_usm_h%ns
    2464                  i = surf_usm_h%i(m)
    2465                  j = surf_usm_h%j(m)
    2466                  k = surf_usm_h%k(m)
    2467                  temp_pf(k,j,i) =  surf_usm_h%frac(0,m)     *                  &
    2468                                    surf_usm_h%emissivity(0,m) +                &
    2469                                    surf_usm_h%frac(1,m)     *                  &
    2470                                    surf_usm_h%emissivity(1,m) +                &
    2471                                    surf_usm_h%frac(2,m)     *                  &
    2472                                    surf_usm_h%emissivity(2,m)
    2473               ENDDO
    2474               DO  l = 0, 3
     2475              IF ( idsint == iup_u )  THEN
     2476                 DO  m = 1, surf_usm_h%ns
     2477                    i = surf_usm_h%i(m)
     2478                    j = surf_usm_h%j(m)
     2479                    k = surf_usm_h%k(m)
     2480                    temp_pf(k,j,i) =  surf_usm_h%frac(0,m)     *                  &
     2481                                      surf_usm_h%emissivity(0,m) +                &
     2482                                      surf_usm_h%frac(1,m)     *                  &
     2483                                      surf_usm_h%emissivity(1,m) +                &
     2484                                      surf_usm_h%frac(2,m)     *                  &
     2485                                      surf_usm_h%emissivity(2,m)
     2486                 ENDDO
     2487              ELSE
     2488                 l = idsidx
    24752489                 DO  m = 1, surf_usm_v(l)%ns
    24762490                    i = surf_usm_v(l)%i(m)
     
    24842498                                     surf_usm_v(l)%emissivity(2,m)
    24852499                 ENDDO
    2486               ENDDO
     2500              ENDIF
    24872501
    24882502          CASE ( 'usm_surfwintrans' )
    24892503!--           transmissivity window tiles
    2490               DO  m = 1, surf_usm_h%ns
    2491                  i = surf_usm_h%i(m)
    2492                  j = surf_usm_h%j(m)
    2493                  k = surf_usm_h%k(m)
    2494                  temp_pf(k,j,i) = surf_usm_h%transmissivity(m)
    2495               ENDDO
    2496               DO  l = 0, 3
     2504              IF ( idsint == iup_u )  THEN
     2505                 DO  m = 1, surf_usm_h%ns
     2506                    i = surf_usm_h%i(m)
     2507                    j = surf_usm_h%j(m)
     2508                    k = surf_usm_h%k(m)
     2509                    temp_pf(k,j,i) = surf_usm_h%transmissivity(m)
     2510                 ENDDO
     2511              ELSE
     2512                 l = idsidx
    24972513                 DO  m = 1, surf_usm_v(l)%ns
    24982514                    i = surf_usm_v(l)%i(m)
     
    25012517                    temp_pf(k,j,i) = surf_usm_v(l)%transmissivity(m)
    25022518                 ENDDO
    2503 
     2519              ENDIF
     2520
     2521          CASE ( 'usm_skyvf' )
     2522!--           sky view factor
     2523              DO isurf = dirstart(ids), dirend(ids)
     2524                 IF ( surfl(id,isurf) == idsint )  THEN
     2525                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = skyvf(isurf)
     2526                 ENDIF
     2527              ENDDO
     2528             
     2529          CASE ( 'usm_skyvft' )
     2530!--           sky view factor
     2531              DO isurf = dirstart(ids), dirend(ids)
     2532                 IF ( surfl(id,isurf) == ids )  THEN
     2533                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = skyvft(isurf)
     2534                 ENDIF
    25042535              ENDDO
    25052536
     
    25182549                             
    25192550                  IF ( surf(ix,isurfs) == is  .AND.  surf(iy,isurfs) == js  .AND.       &
    2520                        surf(iz,isurfs) == ks  .AND.  surf(id,isurfs) == ids )  THEN
     2551                       surf(iz,isurfs) == ks  .AND.  surf(id,isurfs) == idsint )  THEN
    25212552  !--                 correct source surface
    25222553                      temp_pf(surfl(iz,isurflt),surfl(iy,isurflt),surfl(ix,isurflt)) = svf(k,isvf)
     
    25272558!--           array of complete radiation balance
    25282559              IF ( av == 0 )  THEN
    2529                  DO  m = 1, surf_usm_h%ns
    2530                     i = surf_usm_h%i(m)
    2531                     j = surf_usm_h%j(m)
    2532                     k = surf_usm_h%k(m)
    2533                     temp_pf(k,j,i) = surf_usm_h%rad_net_l(m)
    2534                  ENDDO
    2535                  DO  l = 0, 3
     2560                 IF ( idsint == iup_u )  THEN
     2561                    DO  m = 1, surf_usm_h%ns
     2562                       i = surf_usm_h%i(m)
     2563                       j = surf_usm_h%j(m)
     2564                       k = surf_usm_h%k(m)
     2565                       temp_pf(k,j,i) = surf_usm_h%rad_net_l(m)
     2566                    ENDDO
     2567                 ELSE
     2568                    l = idsidx
    25362569                    DO  m = 1, surf_usm_v(l)%ns
    25372570                       i = surf_usm_v(l)%i(m)
     
    25402573                       temp_pf(k,j,i) = surf_usm_v(l)%rad_net_l(m)
    25412574                    ENDDO
    2542                  ENDDO
     2575                 ENDIF
    25432576              ELSE
    2544                  DO  m = 1, surf_usm_h%ns
    2545                     i = surf_usm_h%i(m)
    2546                     j = surf_usm_h%j(m)
    2547                     k = surf_usm_h%k(m)
    2548                     temp_pf(k,j,i) = surf_usm_h%rad_net_av(m)
    2549                  ENDDO
    2550                  DO  l = 0, 3
     2577                 IF ( idsint == iup_u )  THEN
     2578                    DO  m = 1, surf_usm_h%ns
     2579                       i = surf_usm_h%i(m)
     2580                       j = surf_usm_h%j(m)
     2581                       k = surf_usm_h%k(m)
     2582                       temp_pf(k,j,i) = surf_usm_h%rad_net_av(m)
     2583                    ENDDO
     2584                 ELSE
     2585                    l = idsidx
    25512586                    DO  m = 1, surf_usm_v(l)%ns
    25522587                       i = surf_usm_v(l)%i(m)
     
    25552590                       temp_pf(k,j,i) = surf_usm_v(l)%rad_net_av(m)
    25562591                    ENDDO
    2557                  ENDDO
     2592                 ENDIF
    25582593              ENDIF
    25592594
     
    25612596!--           array of sw radiation falling to surface after i-th reflection
    25622597              DO isurf = dirstart(ids), dirend(ids)
    2563                  IF ( surfl(id,isurf) == ids )  THEN
     2598                 IF ( surfl(id,isurf) == idsint )  THEN
    25642599                   IF ( av == 0 )  THEN
    25652600                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinsw(isurf)
     
    25732608!--           array of lw radiation falling to surface after i-th reflection
    25742609              DO isurf = dirstart(ids), dirend(ids)
    2575                  IF ( surfl(id,isurf) == ids )  THEN
     2610                 IF ( surfl(id,isurf) == idsint )  THEN
    25762611                   IF ( av == 0 )  THEN
    25772612                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlw(isurf)
     
    25852620!--           array of direct sw radiation falling to surface from sun
    25862621              DO isurf = dirstart(ids), dirend(ids)
    2587                  IF ( surfl(id,isurf) == ids )  THEN
     2622                 IF ( surfl(id,isurf) == idsint )  THEN
    25882623                   IF ( av == 0 )  THEN
    25892624                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswdir(isurf)
     
    25972632!--           array of difusion sw radiation falling to surface from sky and borders of the domain
    25982633              DO isurf = dirstart(ids), dirend(ids)
    2599                  IF ( surfl(id,isurf) == ids )  THEN
     2634                 IF ( surfl(id,isurf) == idsint )  THEN
    26002635                   IF ( av == 0 )  THEN
    26012636                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswdif(isurf)
     
    26092644!--           array of sw radiation falling to surface from reflections
    26102645              DO isurf = dirstart(ids), dirend(ids)
    2611                  IF ( surfl(id,isurf) == ids )  THEN
     2646                 IF ( surfl(id,isurf) == idsint )  THEN
    26122647                   IF ( av == 0 )  THEN
    26132648                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = &
     
    26192654              ENDDO
    26202655
     2656          CASE ( 'usm_rad_inlwdif' )
     2657!--           array of difusion lw radiation falling to surface from sky and borders of the domain
     2658              DO isurf = dirstart(ids), dirend(ids)
     2659                 IF ( surfl(id,isurf) == idsint )  THEN
     2660                   IF ( av == 0 )  THEN
     2661                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlwdif(isurf)
     2662                   ELSE
     2663                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlwdif_av(isurf)
     2664                   ENDIF
     2665                 ENDIF
     2666              ENDDO
     2667
    26212668          CASE ( 'usm_rad_inlwref' )
    26222669!--           array of lw radiation falling to surface from reflections
    26232670              DO isurf = dirstart(ids), dirend(ids)
    2624                  IF ( surfl(id,isurf) == ids )  THEN
     2671                 IF ( surfl(id,isurf) == idsint )  THEN
    26252672                   IF ( av == 0 )  THEN
    26262673                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlw(isurf) - surfinlwdif(isurf)
     
    26342681!--           array of sw radiation emitted from surface after i-th reflection
    26352682              DO isurf = dirstart(ids), dirend(ids)
    2636                  IF ( surfl(id,isurf) == ids )  THEN
     2683                 IF ( surfl(id,isurf) == idsint )  THEN
    26372684                   IF ( av == 0 )  THEN
    26382685                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfoutsw(isurf)
     
    26462693!--           array of lw radiation emitted from surface after i-th reflection
    26472694              DO isurf = dirstart(ids), dirend(ids)
    2648                  IF ( surfl(id,isurf) == ids )  THEN
     2695                 IF ( surfl(id,isurf) == idsint )  THEN
    26492696                   IF ( av == 0 )  THEN
    26502697                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfoutlw(isurf)
     
    26582705!--           average of array of residua of sw radiation absorbed in surface after last reflection
    26592706              DO isurf = dirstart(ids), dirend(ids)
    2660                  IF ( surfl(id,isurf) == ids )  THEN
     2707                 IF ( surfl(id,isurf) == idsint )  THEN
    26612708                   IF ( av == 0 )  THEN
    26622709                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfins(isurf)
     
    26702717!--           average of array of residua of lw radiation absorbed in surface after last reflection
    26712718              DO isurf = dirstart(ids), dirend(ids)
    2672                  IF ( surfl(id,isurf) == ids )  THEN
     2719                 IF ( surfl(id,isurf) == idsint )  THEN
    26732720                   IF ( av == 0 )  THEN
    26742721                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinl(isurf)
     
    26822729!--           array of heat flux from radiation for surfaces after all reflections
    26832730              IF ( av == 0 )  THEN
    2684                  DO  m = 1, surf_usm_h%ns
    2685                     i = surf_usm_h%i(m)
    2686                     j = surf_usm_h%j(m)
    2687                     k = surf_usm_h%k(m)
    2688                     temp_pf(k,j,i) = surf_usm_h%surfhf(m)
    2689                  ENDDO
    2690                  DO  l = 0, 3
     2731                 IF ( idsint == iup_u )  THEN
     2732                    DO  m = 1, surf_usm_h%ns
     2733                       i = surf_usm_h%i(m)
     2734                       j = surf_usm_h%j(m)
     2735                       k = surf_usm_h%k(m)
     2736                       temp_pf(k,j,i) = surf_usm_h%surfhf(m)
     2737                    ENDDO
     2738                 ELSE
     2739                    l = idsidx
    26912740                    DO  m = 1, surf_usm_v(l)%ns
    26922741                       i = surf_usm_v(l)%i(m)
     
    26952744                       temp_pf(k,j,i) = surf_usm_v(l)%surfhf(m)
    26962745                    ENDDO
    2697                  ENDDO
     2746                 ENDIF
    26982747              ELSE
    2699                  DO  m = 1, surf_usm_h%ns
    2700                     i = surf_usm_h%i(m)
    2701                     j = surf_usm_h%j(m)
    2702                     k = surf_usm_h%k(m)
    2703                     temp_pf(k,j,i) = surf_usm_h%surfhf_av(m)
    2704                  ENDDO
    2705                  DO  l = 0, 3
     2748                 IF ( idsint == iup_u )  THEN
     2749                    DO  m = 1, surf_usm_h%ns
     2750                       i = surf_usm_h%i(m)
     2751                       j = surf_usm_h%j(m)
     2752                       k = surf_usm_h%k(m)
     2753                       temp_pf(k,j,i) = surf_usm_h%surfhf_av(m)
     2754                    ENDDO
     2755                 ELSE
     2756                    l = idsidx
    27062757                    DO  m = 1, surf_usm_v(l)%ns
    27072758                       i = surf_usm_v(l)%i(m)
     
    27102761                       temp_pf(k,j,i) = surf_usm_v(l)%surfhf_av(m)
    27112762                    ENDDO
    2712                  ENDDO
     2763                 ENDIF
    27132764              ENDIF
    27142765 
     
    27162767!--           array of sensible heat flux from surfaces
    27172768              IF ( av == 0 )  THEN
    2718                  DO  m = 1, surf_usm_h%ns
    2719                     i = surf_usm_h%i(m)
    2720                     j = surf_usm_h%j(m)
    2721                     k = surf_usm_h%k(m)
    2722                     temp_pf(k,j,i) = surf_usm_h%wshf_eb(m)
    2723                  ENDDO
    2724                  DO  l = 0, 3
     2769                 IF ( idsint == iup_u )  THEN
     2770                    DO  m = 1, surf_usm_h%ns
     2771                       i = surf_usm_h%i(m)
     2772                       j = surf_usm_h%j(m)
     2773                       k = surf_usm_h%k(m)
     2774                       temp_pf(k,j,i) = surf_usm_h%wshf_eb(m)
     2775                    ENDDO
     2776                 ELSE
     2777                    l = idsidx
    27252778                    DO  m = 1, surf_usm_v(l)%ns
    27262779                       i = surf_usm_v(l)%i(m)
     
    27292782                       temp_pf(k,j,i) = surf_usm_v(l)%wshf_eb(m)
    27302783                    ENDDO
    2731                  ENDDO
     2784                 ENDIF
    27322785              ELSE
    2733                  DO  m = 1, surf_usm_h%ns
    2734                     i = surf_usm_h%i(m)
    2735                     j = surf_usm_h%j(m)
    2736                     k = surf_usm_h%k(m)
    2737                     temp_pf(k,j,i) = surf_usm_h%wshf_eb_av(m)
    2738                  ENDDO
    2739                  DO  l = 0, 3
     2786                 IF ( idsint == iup_u )  THEN
     2787                    DO  m = 1, surf_usm_h%ns
     2788                       i = surf_usm_h%i(m)
     2789                       j = surf_usm_h%j(m)
     2790                       k = surf_usm_h%k(m)
     2791                       temp_pf(k,j,i) = surf_usm_h%wshf_eb_av(m)
     2792                    ENDDO
     2793                 ELSE
     2794                    l = idsidx
    27402795                    DO  m = 1, surf_usm_v(l)%ns
    27412796                       i = surf_usm_v(l)%i(m)
     
    27442799                       temp_pf(k,j,i) = surf_usm_v(l)%wshf_eb_av(m)
    27452800                    ENDDO
    2746                  ENDDO
     2801                 ENDIF
    27472802              ENDIF
    27482803
     
    27512806!--           array of heat flux from ground (land, wall, roof)
    27522807              IF ( av == 0 )  THEN
    2753                  DO  m = 1, surf_usm_h%ns
    2754                     i = surf_usm_h%i(m)
    2755                     j = surf_usm_h%j(m)
    2756                     k = surf_usm_h%k(m)
    2757                     temp_pf(k,j,i) = surf_usm_h%wghf_eb(m)
    2758                  ENDDO
    2759                  DO  l = 0, 3
     2808                 IF ( idsint == iup_u )  THEN
     2809                    DO  m = 1, surf_usm_h%ns
     2810                       i = surf_usm_h%i(m)
     2811                       j = surf_usm_h%j(m)
     2812                       k = surf_usm_h%k(m)
     2813                       temp_pf(k,j,i) = surf_usm_h%wghf_eb(m)
     2814                    ENDDO
     2815                 ELSE
     2816                    l = idsidx
    27602817                    DO  m = 1, surf_usm_v(l)%ns
    27612818                       i = surf_usm_v(l)%i(m)
     
    27642821                       temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb(m)
    27652822                    ENDDO
    2766                  ENDDO
     2823                 ENDIF
    27672824              ELSE
    2768                  DO  m = 1, surf_usm_h%ns
    2769                     i = surf_usm_h%i(m)
    2770                     j = surf_usm_h%j(m)
    2771                     k = surf_usm_h%k(m)
    2772                     temp_pf(k,j,i) = surf_usm_h%wghf_eb_av(m)
    2773                  ENDDO
    2774                  DO  l = 0, 3
     2825                 IF ( idsint == iup_u )  THEN
     2826                    DO  m = 1, surf_usm_h%ns
     2827                       i = surf_usm_h%i(m)
     2828                       j = surf_usm_h%j(m)
     2829                       k = surf_usm_h%k(m)
     2830                       temp_pf(k,j,i) = surf_usm_h%wghf_eb_av(m)
     2831                    ENDDO
     2832                 ELSE
     2833                    l = idsidx
    27752834                    DO  m = 1, surf_usm_v(l)%ns
    27762835                       i = surf_usm_v(l)%i(m)
     
    27792838                       temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_av(m)
    27802839                    ENDDO
    2781                  ENDDO
     2840                 ENDIF
    27822841              ENDIF
    27832842
     
    27862845
    27872846              IF ( av == 0 )  THEN
    2788                  DO  m = 1, surf_usm_h%ns
    2789                     i = surf_usm_h%i(m)
    2790                     j = surf_usm_h%j(m)
    2791                     k = surf_usm_h%k(m)
    2792                     temp_pf(k,j,i) = surf_usm_h%wghf_eb_window(m)
    2793                  ENDDO
    2794                  DO  l = 0, 3
     2847                 IF ( idsint == iup_u )  THEN
     2848                    DO  m = 1, surf_usm_h%ns
     2849                       i = surf_usm_h%i(m)
     2850                       j = surf_usm_h%j(m)
     2851                       k = surf_usm_h%k(m)
     2852                       temp_pf(k,j,i) = surf_usm_h%wghf_eb_window(m)
     2853                    ENDDO
     2854                 ELSE
     2855                    l = idsidx
    27952856                    DO  m = 1, surf_usm_v(l)%ns
    27962857                       i = surf_usm_v(l)%i(m)
     
    27992860                       temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_window(m)
    28002861                    ENDDO
    2801                  ENDDO
     2862                 ENDIF
    28022863              ELSE
    2803                  DO  m = 1, surf_usm_h%ns
    2804                     i = surf_usm_h%i(m)
    2805                     j = surf_usm_h%j(m)
    2806                     k = surf_usm_h%k(m)
    2807                     temp_pf(k,j,i) = surf_usm_h%wghf_eb_window_av(m)
    2808                  ENDDO
    2809                  DO  l = 0, 3
     2864                 IF ( idsint == iup_u )  THEN
     2865                    DO  m = 1, surf_usm_h%ns
     2866                       i = surf_usm_h%i(m)
     2867                       j = surf_usm_h%j(m)
     2868                       k = surf_usm_h%k(m)
     2869                       temp_pf(k,j,i) = surf_usm_h%wghf_eb_window_av(m)
     2870                    ENDDO
     2871                 ELSE
     2872                    l = idsidx
    28102873                    DO  m = 1, surf_usm_v(l)%ns
    28112874                       i = surf_usm_v(l)%i(m)
     
    28142877                       temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_window_av(m)
    28152878                    ENDDO
    2816                  ENDDO
     2879                 ENDIF
    28172880              ENDIF
    28182881
     
    28212884
    28222885              IF ( av == 0 )  THEN
    2823                  DO  m = 1, surf_usm_h%ns
    2824                     i = surf_usm_h%i(m)
    2825                     j = surf_usm_h%j(m)
    2826                     k = surf_usm_h%k(m)
    2827                     temp_pf(k,j,i) = surf_usm_h%wghf_eb_green(m)
    2828                  ENDDO
    2829                  DO  l = 0, 3
     2886                 IF ( idsint == iup_u )  THEN
     2887                    DO  m = 1, surf_usm_h%ns
     2888                       i = surf_usm_h%i(m)
     2889                       j = surf_usm_h%j(m)
     2890                       k = surf_usm_h%k(m)
     2891                       temp_pf(k,j,i) = surf_usm_h%wghf_eb_green(m)
     2892                    ENDDO
     2893                 ELSE
     2894                    l = idsidx
    28302895                    DO  m = 1, surf_usm_v(l)%ns
    28312896                       i = surf_usm_v(l)%i(m)
     
    28342899                       temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_green(m)
    28352900                    ENDDO
    2836                  ENDDO
     2901                 ENDIF
    28372902              ELSE
    2838                  DO  m = 1, surf_usm_h%ns
    2839                     i = surf_usm_h%i(m)
    2840                     j = surf_usm_h%j(m)
    2841                     k = surf_usm_h%k(m)
    2842                     temp_pf(k,j,i) = surf_usm_h%wghf_eb_green_av(m)
    2843                  ENDDO
    2844                  DO  l = 0, 3
     2903                 IF ( idsint == iup_u )  THEN
     2904                    DO  m = 1, surf_usm_h%ns
     2905                       i = surf_usm_h%i(m)
     2906                       j = surf_usm_h%j(m)
     2907                       k = surf_usm_h%k(m)
     2908                       temp_pf(k,j,i) = surf_usm_h%wghf_eb_green_av(m)
     2909                    ENDDO
     2910                 ELSE
     2911                    l = idsidx
    28452912                    DO  m = 1, surf_usm_v(l)%ns
    28462913                       i = surf_usm_v(l)%i(m)
     
    28492916                       temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_green_av(m)
    28502917                    ENDDO
    2851                  ENDDO
     2918                 ENDIF
    28522919              ENDIF
    28532920
     
    28552922!--           array of heat flux from indoor ground (land, wall, roof)
    28562923              IF ( av == 0 )  THEN
    2857                  DO  m = 1, surf_usm_h%ns
    2858                     i = surf_usm_h%i(m)
    2859                     j = surf_usm_h%j(m)
    2860                     k = surf_usm_h%k(m)
    2861                     temp_pf(k,j,i) = surf_usm_h%iwghf_eb(m)
    2862                  ENDDO
    2863                  DO  l = 0, 3
     2924                 IF ( idsint == iup_u )  THEN
     2925                    DO  m = 1, surf_usm_h%ns
     2926                       i = surf_usm_h%i(m)
     2927                       j = surf_usm_h%j(m)
     2928                       k = surf_usm_h%k(m)
     2929                       temp_pf(k,j,i) = surf_usm_h%iwghf_eb(m)
     2930                    ENDDO
     2931                 ELSE
     2932                    l = idsidx
    28642933                    DO  m = 1, surf_usm_v(l)%ns
    28652934                       i = surf_usm_v(l)%i(m)
     
    28682937                       temp_pf(k,j,i) = surf_usm_v(l)%iwghf_eb(m)
    28692938                    ENDDO
    2870                  ENDDO
     2939                 ENDIF
    28712940              ELSE
    2872                  DO  m = 1, surf_usm_h%ns
    2873                     i = surf_usm_h%i(m)
    2874                     j = surf_usm_h%j(m)
    2875                     k = surf_usm_h%k(m)
    2876                     temp_pf(k,j,i) = surf_usm_h%iwghf_eb_av(m)
    2877                  ENDDO
    2878                  DO  l = 0, 3
     2941                 IF ( idsint == iup_u )  THEN
     2942                    DO  m = 1, surf_usm_h%ns
     2943                       i = surf_usm_h%i(m)
     2944                       j = surf_usm_h%j(m)
     2945                       k = surf_usm_h%k(m)
     2946                       temp_pf(k,j,i) = surf_usm_h%iwghf_eb_av(m)
     2947                    ENDDO
     2948                 ELSE
     2949                    l = idsidx
    28792950                    DO  m = 1, surf_usm_v(l)%ns
    28802951                       i = surf_usm_v(l)%i(m)
     
    28832954                       temp_pf(k,j,i) = surf_usm_v(l)%iwghf_eb_av(m)
    28842955                    ENDDO
    2885                  ENDDO
     2956                 ENDIF
    28862957              ENDIF
    28872958
     
    28902961
    28912962              IF ( av == 0 )  THEN
    2892                  DO  m = 1, surf_usm_h%ns
    2893                     i = surf_usm_h%i(m)
    2894                     j = surf_usm_h%j(m)
    2895                     k = surf_usm_h%k(m)
    2896                     temp_pf(k,j,i) = surf_usm_h%iwghf_eb_window(m)
    2897                  ENDDO
    2898                  DO  l = 0, 3
     2963                 IF ( idsint == iup_u )  THEN
     2964                    DO  m = 1, surf_usm_h%ns
     2965                       i = surf_usm_h%i(m)
     2966                       j = surf_usm_h%j(m)
     2967                       k = surf_usm_h%k(m)
     2968                       temp_pf(k,j,i) = surf_usm_h%iwghf_eb_window(m)
     2969                    ENDDO
     2970                 ELSE
     2971                    l = idsidx
    28992972                    DO  m = 1, surf_usm_v(l)%ns
    29002973                       i = surf_usm_v(l)%i(m)
     
    29032976                       temp_pf(k,j,i) = surf_usm_v(l)%iwghf_eb_window(m)
    29042977                    ENDDO
    2905                  ENDDO
     2978                 ENDIF
    29062979              ELSE
    2907                  DO  m = 1, surf_usm_h%ns
    2908                     i = surf_usm_h%i(m)
    2909                     j = surf_usm_h%j(m)
    2910                     k = surf_usm_h%k(m)
    2911                     temp_pf(k,j,i) = surf_usm_h%iwghf_eb_window_av(m)
    2912                  ENDDO
    2913                  DO  l = 0, 3
     2980                 IF ( idsint == iup_u )  THEN
     2981                    DO  m = 1, surf_usm_h%ns
     2982                       i = surf_usm_h%i(m)
     2983                       j = surf_usm_h%j(m)
     2984                       k = surf_usm_h%k(m)
     2985                       temp_pf(k,j,i) = surf_usm_h%iwghf_eb_window_av(m)
     2986                    ENDDO
     2987                 ELSE
     2988                    l = idsidx
    29142989                    DO  m = 1, surf_usm_v(l)%ns
    29152990                       i = surf_usm_v(l)%i(m)
     
    29182993                       temp_pf(k,j,i) = surf_usm_v(l)%iwghf_eb_window_av(m)
    29192994                    ENDDO
    2920                  ENDDO
     2995                 ENDIF
    29212996              ENDIF
    29222997             
     
    29242999!--           surface temperature for surfaces
    29253000              IF ( av == 0 )  THEN
    2926                  DO  m = 1, surf_usm_h%ns
    2927                     i = surf_usm_h%i(m)
    2928                     j = surf_usm_h%j(m)
    2929                     k = surf_usm_h%k(m)
    2930                     temp_pf(k,j,i) = t_surf_h(m)
    2931                  ENDDO
    2932                  DO  l = 0, 3
     3001                 IF ( idsint == iup_u )  THEN
     3002                    DO  m = 1, surf_usm_h%ns
     3003                       i = surf_usm_h%i(m)
     3004                       j = surf_usm_h%j(m)
     3005                       k = surf_usm_h%k(m)
     3006                       temp_pf(k,j,i) = t_surf_h(m)
     3007                    ENDDO
     3008                 ELSE
     3009                    l = idsidx
    29333010                    DO  m = 1, surf_usm_v(l)%ns
    29343011                       i = surf_usm_v(l)%i(m)
     
    29373014                       temp_pf(k,j,i) = t_surf_v(l)%t(m)
    29383015                    ENDDO
    2939                  ENDDO
     3016                 ENDIF
    29403017              ELSE
    2941                  DO  m = 1, surf_usm_h%ns
    2942                     i = surf_usm_h%i(m)
    2943                     j = surf_usm_h%j(m)
    2944                     k = surf_usm_h%k(m)
    2945                     temp_pf(k,j,i) = surf_usm_h%t_surf_av(m)
    2946                  ENDDO
    2947                  DO  l = 0, 3
     3018                 IF ( idsint == iup_u )  THEN
     3019                    DO  m = 1, surf_usm_h%ns
     3020                       i = surf_usm_h%i(m)
     3021                       j = surf_usm_h%j(m)
     3022                       k = surf_usm_h%k(m)
     3023                       temp_pf(k,j,i) = surf_usm_h%t_surf_av(m)
     3024                    ENDDO
     3025                 ELSE
     3026                    l = idsidx
    29483027                    DO  m = 1, surf_usm_v(l)%ns
    29493028                       i = surf_usm_v(l)%i(m)
     
    29523031                       temp_pf(k,j,i) = surf_usm_v(l)%t_surf_av(m)
    29533032                    ENDDO
    2954                  ENDDO
     3033                 ENDIF
    29553034              ENDIF
    29563035             
     
    29593038
    29603039              IF ( av == 0 )  THEN
    2961                  DO  m = 1, surf_usm_h%ns
    2962                     i = surf_usm_h%i(m)
    2963                     j = surf_usm_h%j(m)
    2964                     k = surf_usm_h%k(m)
    2965                     temp_pf(k,j,i) = t_surf_window_h(m)
    2966                  ENDDO
    2967                  DO  l = 0, 3
     3040                 IF ( idsint == iup_u )  THEN
     3041                    DO  m = 1, surf_usm_h%ns
     3042                       i = surf_usm_h%i(m)
     3043                       j = surf_usm_h%j(m)
     3044                       k = surf_usm_h%k(m)
     3045                       temp_pf(k,j,i) = t_surf_window_h(m)
     3046                    ENDDO
     3047                 ELSE
     3048                    l = idsidx
    29683049                    DO  m = 1, surf_usm_v(l)%ns
    29693050                       i = surf_usm_v(l)%i(m)
     
    29723053                       temp_pf(k,j,i) = t_surf_window_v(l)%t(m)
    29733054                    ENDDO
    2974                  ENDDO
     3055                 ENDIF
    29753056
    29763057              ELSE
    2977                  DO  m = 1, surf_usm_h%ns
    2978                     i = surf_usm_h%i(m)
    2979                     j = surf_usm_h%j(m)
    2980                     k = surf_usm_h%k(m)
    2981                     temp_pf(k,j,i) = surf_usm_h%t_surf_window_av(m)
    2982                  ENDDO
    2983                  DO  l = 0, 3
     3058                 IF ( idsint == iup_u )  THEN
     3059                    DO  m = 1, surf_usm_h%ns
     3060                       i = surf_usm_h%i(m)
     3061                       j = surf_usm_h%j(m)
     3062                       k = surf_usm_h%k(m)
     3063                       temp_pf(k,j,i) = surf_usm_h%t_surf_window_av(m)
     3064                    ENDDO
     3065                 ELSE
     3066                    l = idsidx
    29843067                    DO  m = 1, surf_usm_v(l)%ns
    29853068                       i = surf_usm_v(l)%i(m)
     
    29893072                    ENDDO
    29903073
    2991                  ENDDO
     3074                 ENDIF
    29923075
    29933076              ENDIF
     
    29973080
    29983081              IF ( av == 0 )  THEN
    2999                  DO  m = 1, surf_usm_h%ns
    3000                     i = surf_usm_h%i(m)
    3001                     j = surf_usm_h%j(m)
    3002                     k = surf_usm_h%k(m)
    3003                     temp_pf(k,j,i) = t_surf_green_h(m)
    3004                  ENDDO
    3005                  DO  l = 0, 3
     3082                 IF ( idsint == iup_u )  THEN
     3083                    DO  m = 1, surf_usm_h%ns
     3084                       i = surf_usm_h%i(m)
     3085                       j = surf_usm_h%j(m)
     3086                       k = surf_usm_h%k(m)
     3087                       temp_pf(k,j,i) = t_surf_green_h(m)
     3088                    ENDDO
     3089                 ELSE
     3090                    l = idsidx
    30063091                    DO  m = 1, surf_usm_v(l)%ns
    30073092                       i = surf_usm_v(l)%i(m)
     
    30103095                       temp_pf(k,j,i) = t_surf_green_v(l)%t(m)
    30113096                    ENDDO
    3012                  ENDDO
     3097                 ENDIF
    30133098
    30143099              ELSE
    3015                  DO  m = 1, surf_usm_h%ns
    3016                     i = surf_usm_h%i(m)
    3017                     j = surf_usm_h%j(m)
    3018                     k = surf_usm_h%k(m)
    3019                     temp_pf(k,j,i) = surf_usm_h%t_surf_green_av(m)
    3020                  ENDDO
    3021                  DO  l = 0, 3
     3100                 IF ( idsint == iup_u )  THEN
     3101                    DO  m = 1, surf_usm_h%ns
     3102                       i = surf_usm_h%i(m)
     3103                       j = surf_usm_h%j(m)
     3104                       k = surf_usm_h%k(m)
     3105                       temp_pf(k,j,i) = surf_usm_h%t_surf_green_av(m)
     3106                    ENDDO
     3107                 ELSE
     3108                    l = idsidx
    30223109                    DO  m = 1, surf_usm_v(l)%ns
    30233110                       i = surf_usm_v(l)%i(m)
     
    30273114                    ENDDO
    30283115
    3029                  ENDDO
     3116                 ENDIF
    30303117
    30313118              ENDIF
     
    30353122
    30363123              IF ( av == 0 )  THEN
    3037                  DO  m = 1, surf_usm_h%ns
    3038                     i = surf_usm_h%i(m)
    3039                     j = surf_usm_h%j(m)
    3040                     k = surf_usm_h%k(m)
    3041                     temp_pf(k,j,i) = t_surf_10cm_h(m)
    3042                  ENDDO
    3043                  DO  l = 0, 3
     3124                 IF ( idsint == iup_u )  THEN
     3125                    DO  m = 1, surf_usm_h%ns
     3126                       i = surf_usm_h%i(m)
     3127                       j = surf_usm_h%j(m)
     3128                       k = surf_usm_h%k(m)
     3129                       temp_pf(k,j,i) = t_surf_10cm_h(m)
     3130                    ENDDO
     3131                 ELSE
     3132                    l = idsidx
    30443133                    DO  m = 1, surf_usm_v(l)%ns
    30453134                       i = surf_usm_v(l)%i(m)
     
    30483137                       temp_pf(k,j,i) = t_surf_10cm_v(l)%t(m)
    30493138                    ENDDO
    3050                  ENDDO
     3139                 ENDIF
    30513140
    30523141              ELSE
    3053                  DO  m = 1, surf_usm_h%ns
    3054                     i = surf_usm_h%i(m)
    3055                     j = surf_usm_h%j(m)
    3056                     k = surf_usm_h%k(m)
    3057                     temp_pf(k,j,i) = surf_usm_h%t_surf_10cm_av(m)
    3058                  ENDDO
    3059                  DO  l = 0, 3
     3142                 IF ( idsint == iup_u )  THEN
     3143                    DO  m = 1, surf_usm_h%ns
     3144                       i = surf_usm_h%i(m)
     3145                       j = surf_usm_h%j(m)
     3146                       k = surf_usm_h%k(m)
     3147                       temp_pf(k,j,i) = surf_usm_h%t_surf_10cm_av(m)
     3148                    ENDDO
     3149                 ELSE
     3150                    l = idsidx
    30603151                    DO  m = 1, surf_usm_v(l)%ns
    30613152                       i = surf_usm_v(l)%i(m)
     
    30653156                    ENDDO
    30663157
    3067                  ENDDO
     3158                 ENDIF
    30683159
    30693160              ENDIF
     
    30733164!--           wall temperature for  iwl layer of walls and land
    30743165              IF ( av == 0 )  THEN
    3075                  DO  m = 1, surf_usm_h%ns
    3076                     i = surf_usm_h%i(m)
    3077                     j = surf_usm_h%j(m)
    3078                     k = surf_usm_h%k(m)
    3079                     temp_pf(k,j,i) = t_wall_h(iwl,m)
    3080                  ENDDO
    3081                  DO  l = 0, 3
     3166                 IF ( idsint == iup_u )  THEN
     3167                    DO  m = 1, surf_usm_h%ns
     3168                       i = surf_usm_h%i(m)
     3169                       j = surf_usm_h%j(m)
     3170                       k = surf_usm_h%k(m)
     3171                       temp_pf(k,j,i) = t_wall_h(iwl,m)
     3172                    ENDDO
     3173                 ELSE
     3174                    l = idsidx
    30823175                    DO  m = 1, surf_usm_v(l)%ns
    30833176                       i = surf_usm_v(l)%i(m)
     
    30863179                       temp_pf(k,j,i) = t_wall_v(l)%t(iwl,m)
    30873180                    ENDDO
    3088                  ENDDO
     3181                 ENDIF
    30893182              ELSE
    3090                  DO  m = 1, surf_usm_h%ns
    3091                     i = surf_usm_h%i(m)
    3092                     j = surf_usm_h%j(m)
    3093                     k = surf_usm_h%k(m)
    3094                     temp_pf(k,j,i) = surf_usm_h%t_wall_av(iwl,m)
    3095                  ENDDO
    3096                  DO  l = 0, 3
     3183                 IF ( idsint == iup_u )  THEN
     3184                    DO  m = 1, surf_usm_h%ns
     3185                       i = surf_usm_h%i(m)
     3186                       j = surf_usm_h%j(m)
     3187                       k = surf_usm_h%k(m)
     3188                       temp_pf(k,j,i) = surf_usm_h%t_wall_av(iwl,m)
     3189                    ENDDO
     3190                 ELSE
     3191                    l = idsidx
    30973192                    DO  m = 1, surf_usm_v(l)%ns
    30983193                       i = surf_usm_v(l)%i(m)
     
    31013196                       temp_pf(k,j,i) = surf_usm_v(l)%t_wall_av(iwl,m)
    31023197                    ENDDO
    3103                  ENDDO
     3198                 ENDIF
    31043199              ENDIF
    31053200             
     
    31073202!--           window temperature for iwl layer of walls and land
    31083203              IF ( av == 0 )  THEN
    3109                  DO  m = 1, surf_usm_h%ns
    3110                     i = surf_usm_h%i(m)
    3111                     j = surf_usm_h%j(m)
    3112                     k = surf_usm_h%k(m)
    3113                     temp_pf(k,j,i) = t_window_h(iwl,m)
    3114                  ENDDO
    3115                  DO  l = 0, 3
     3204                 IF ( idsint == iup_u )  THEN
     3205                    DO  m = 1, surf_usm_h%ns
     3206                       i = surf_usm_h%i(m)
     3207                       j = surf_usm_h%j(m)
     3208                       k = surf_usm_h%k(m)
     3209                       temp_pf(k,j,i) = t_window_h(iwl,m)
     3210                    ENDDO
     3211                 ELSE
     3212                    l = idsidx
    31163213                    DO  m = 1, surf_usm_v(l)%ns
    31173214                       i = surf_usm_v(l)%i(m)
     
    31203217                       temp_pf(k,j,i) = t_window_v(l)%t(iwl,m)
    31213218                    ENDDO
    3122                  ENDDO
     3219                 ENDIF
    31233220              ELSE
    3124                  DO  m = 1, surf_usm_h%ns
    3125                     i = surf_usm_h%i(m)
    3126                     j = surf_usm_h%j(m)
    3127                     k = surf_usm_h%k(m)
    3128                     temp_pf(k,j,i) = surf_usm_h%t_window_av(iwl,m)
    3129                  ENDDO
    3130                  DO  l = 0, 3
     3221                 IF ( idsint == iup_u )  THEN
     3222                    DO  m = 1, surf_usm_h%ns
     3223                       i = surf_usm_h%i(m)
     3224                       j = surf_usm_h%j(m)
     3225                       k = surf_usm_h%k(m)
     3226                       temp_pf(k,j,i) = surf_usm_h%t_window_av(iwl,m)
     3227                    ENDDO
     3228                 ELSE
     3229                    l = idsidx
    31313230                    DO  m = 1, surf_usm_v(l)%ns
    31323231                       i = surf_usm_v(l)%i(m)
     
    31353234                       temp_pf(k,j,i) = surf_usm_v(l)%t_window_av(iwl,m)
    31363235                    ENDDO
    3137                  ENDDO
     3236                 ENDIF
    31383237              ENDIF
    31393238
     
    31413240!--           green temperature for  iwl layer of walls and land
    31423241              IF ( av == 0 )  THEN
    3143                  DO  m = 1, surf_usm_h%ns
    3144                     i = surf_usm_h%i(m)
    3145                     j = surf_usm_h%j(m)
    3146                     k = surf_usm_h%k(m)
    3147                     temp_pf(k,j,i) = t_green_h(iwl,m)
    3148                  ENDDO
    3149                  DO  l = 0, 3
     3242                 IF ( idsint == iup_u )  THEN
     3243                    DO  m = 1, surf_usm_h%ns
     3244                       i = surf_usm_h%i(m)
     3245                       j = surf_usm_h%j(m)
     3246                       k = surf_usm_h%k(m)
     3247                       temp_pf(k,j,i) = t_green_h(iwl,m)
     3248                    ENDDO
     3249                 ELSE
     3250                    l = idsidx
    31503251                    DO  m = 1, surf_usm_v(l)%ns
    31513252                       i = surf_usm_v(l)%i(m)
     
    31543255                       temp_pf(k,j,i) = t_green_v(l)%t(iwl,m)
    31553256                    ENDDO
    3156                  ENDDO
     3257                 ENDIF
    31573258              ELSE
    3158                  DO  m = 1, surf_usm_h%ns
    3159                     i = surf_usm_h%i(m)
    3160                     j = surf_usm_h%j(m)
    3161                     k = surf_usm_h%k(m)
    3162                     temp_pf(k,j,i) = surf_usm_h%t_green_av(iwl,m)
    3163                  ENDDO
    3164                  DO  l = 0, 3
     3259                 IF ( idsint == iup_u )  THEN
     3260                    DO  m = 1, surf_usm_h%ns
     3261                       i = surf_usm_h%i(m)
     3262                       j = surf_usm_h%j(m)
     3263                       k = surf_usm_h%k(m)
     3264                       temp_pf(k,j,i) = surf_usm_h%t_green_av(iwl,m)
     3265                    ENDDO
     3266                 ELSE
     3267                    l = idsidx
    31653268                    DO  m = 1, surf_usm_v(l)%ns
    31663269                       i = surf_usm_v(l)%i(m)
     
    31693272                       temp_pf(k,j,i) = surf_usm_v(l)%t_green_av(iwl,m)
    31703273                    ENDDO
    3171                  ENDDO
     3274                 ENDIF
    31723275              ENDIF
    31733276
     
    31803283!
    31813284!--     Rearrange dimensions for NetCDF output
     3285!--     FIXME: this may generate FPE overflow upon conversion from DP to SP
    31823286        DO  j = nys, nyn
    31833287            DO  i = nxl, nxr
    31843288                DO  k = nzb_do, nzt_do
     3289!                   print*, j,nys,nyn,i,nxl,nxr,k,nzb_do,nzt_do,local_pf(i,j,k),temp_pf(k,j,i)
    31853290                    local_pf(i,j,k) = temp_pf(k,j,i)
    31863291                ENDDO
     
    32283333             var(1:7) == 'usm_dif'  .OR.  var(1:11) == 'usm_surfcat'  .OR.                  &
    32293334             var(1:11) == 'usm_surfalb'  .OR.  var(1:12) == 'usm_surfemis'  .OR.            &
    3230              var(1:16) == 'usm_surfwintrans' )  THEN
     3335             var(1:16) == 'usm_surfwintrans'  .OR.                                          &
     3336             var(1:9) == 'usm_skyvf' .OR. var(1:9) == 'usm_skyvft' ) THEN
    32313337
    32323338            found = .TRUE.
     
    34163522        INTEGER(iwp) ::  st                  !< dummy 
    34173523
    3418         REAL(wp)     ::  c, d, tin, twin, exn
     3524        REAL(wp)     ::  c, d, tin, twin
    34193525        REAL(wp)     ::  ground_floor_level_l         !< local height of ground floor level
    34203526        REAL(wp)     ::  z_agl                        !< height above ground
     3527        REAL(wp), DIMENSION(nzb:nzt)   ::  exn        !< value of the Exner function in layers
    34213528
    34223529!
     
    44024509            CALL usm_read_anthropogenic_heat()
    44034510        ENDIF
    4404                
     4511
    44054512        IF ( plant_canopy )  THEN
    44064513           
     
    44194526       
    44204527!--         Calculate initial surface temperature from pt of adjacent gridbox
    4421             exn = ( surface_pressure / 1000.0_wp )**0.286_wp
     4528#if ! defined( __nopointer )
     4529            exn(nzb:nzt) = (hyp(nzb:nzt) / 100000.0_wp )**0.286_wp          !< Exner function
     4530#endif
    44224531
    44234532!
     
    44304539               k = surf_usm_h%k(m)
    44314540
    4432                t_surf_h(m) = pt(k,j,i) * exn
    4433                t_surf_window_h(m) = pt(k,j,i) * exn
    4434                t_surf_green_h(m) = pt(k,j,i) * exn
    4435                surf_usm_h%pt_surface(m) = pt(k,j,i) * exn
     4541               t_surf_h(m) = pt(k,j,i) * exn(k)
     4542               t_surf_window_h(m) = pt(k,j,i) * exn(k)
     4543               t_surf_green_h(m) = pt(k,j,i) * exn(k)
     4544               surf_usm_h%pt_surface(m) = pt(k,j,i) * exn(k)
    44364545            ENDDO
    44374546!
     
    44434552                  k = surf_usm_v(l)%k(m)
    44444553
    4445                   t_surf_v(l)%t(m) = pt(k,j,i) * exn
    4446                   t_surf_window_v(l)%t(m) = pt(k,j,i) * exn
    4447                   t_surf_green_v(l)%t(m) = pt(k,j,i) * exn
    4448                   surf_usm_v(l)%pt_surface(m) = pt(k,j,i) * exn
     4554                  t_surf_v(l)%t(m) = pt(k,j,i) * exn(k)
     4555                  t_surf_window_v(l)%t(m) = pt(k,j,i) * exn(k)
     4556                  t_surf_green_v(l)%t(m) = pt(k,j,i) * exn(k)
     4557                  surf_usm_v(l)%pt_surface(m) = pt(k,j,i) * exn(k)
    44494558               ENDDO
    44504559            ENDDO
     
    44964605               ENDDO
    44974606            ENDDO
    4498 
     4607        ELSE
     4608!--         If specified, replace constant wall temperatures with fully 3D values from file
     4609            IF ( read_wall_temp_3d )  CALL usm_read_wall_temperature()
     4610!
    44994611        ENDIF
    45004612       
     
    45194631        t_green_h_p = t_green_h
    45204632        t_green_v_p = t_green_v
     4633
     4634!--     Adjust radiative fluxes for urban surface at model start
     4635        !CALL radiation_interaction
     4636!--     TODO: interaction should be called once before first output,
     4637!--     that is not yet possible.
    45214638       
    45224639        CALL cpu_log( log_point_s(78), 'usm_init', 'stop' )
     
    50295146                           building_type,                                      &
    50305147                           land_category,                                      &
    5031                            pedestrant_category,                                &
    5032                            ra_horiz_coef,                                      &
     5148                           naheatlayers,                                       &
     5149                           pedestrian_category,                                &
     5150                           roughness_concrete,                                 &
     5151                           read_wall_temp_3d,                                  &
    50335152                           roof_category,                                      &
    50345153                           urban_surface,                                      &
    50355154                           usm_anthropogenic_heat,                             &
    50365155                           usm_material_model,                                 &
    5037                            usm_lad_rma,                                        &
    50385156                           wall_category,                                      &
    5039                            indoor_model
    5040 
    5041        line = ' '
     5157                           indoor_model,                                       &
     5158                           wall_inner_temperature,                             &
     5159                           roof_inner_temperature,                             &
     5160                           soil_inner_temperature,                             &
     5161                           window_inner_temperature
    50425162
    50435163!
     
    51345254    SUBROUTINE usm_read_anthropogenic_heat
    51355255   
    5136         INTEGER(iwp)                  :: i,j,ii
     5256        INTEGER(iwp)                  :: i,j,k,ii
    51375257        REAL(wp)                      :: heat
    5138        
     5258
    51395259!--     allocation of array of sources of anthropogenic heat and their diural profile
    5140         ALLOCATE( aheat(nys:nyn,nxl:nxr) )
    5141         ALLOCATE( aheatprof(0:24) )
     5260        ALLOCATE( aheat(naheatlayers,nys:nyn,nxl:nxr) )
     5261        ALLOCATE( aheatprof(naheatlayers,0:24) )
    51425262
    51435263!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     
    51545274                j = 0
    51555275                DO
    5156                     READ( 151, *, err=12, end=13 )  i, j, heat
     5276                    READ( 151, *, err=12, end=13 )  i, j, k, heat
    51575277                    IF ( i >= nxl  .AND.  i <= nxr  .AND.  j >= nys  .AND.  j <= nyn )  THEN
    5158 !--                     write heat into the array
    5159                         aheat(j,i) = heat
     5278                        IF ( k <= naheatlayers  .AND.  k > get_topography_top_index_ji( j, i, 's' ) )  THEN
     5279!--                         write heat into the array
     5280                            aheat(k,j,i) = heat
     5281                        ENDIF
    51605282                    ENDIF
    51615283                    CYCLE
     
    51865308                i = 0
    51875309                DO
    5188                     READ( 151, *, err=22, end=23 )  i, heat
    5189                     IF ( i >= 0  .AND.  i <= 24 )  THEN
     5310                    READ( 151, *, err=22, end=23 )  i, k, heat
     5311                    IF ( i >= 0  .AND.  i <= 24  .AND.  k <= naheatlayers )  THEN
    51905312!--                     write heat into the array
    5191                         aheatprof(i) = heat
     5313                        aheatprof(k,i) = heat
    51925314                    ENDIF
    51935315                    CYCLE
     
    51965318                    CALL message( 'usm_read_anthropogenic_heat', 'PA0517', 0, 1, 0, 6, 0 )
    51975319                ENDDO
    5198                 aheatprof(24) = aheatprof(0)
     5320                aheatprof(:,24) = aheatprof(:,0)
    51995321 23             CLOSE(151)
    52005322                CYCLE
     
    52125334
    52135335!------------------------------------------------------------------------------!
    5214 !
    52155336! Description:
    52165337! ------------
    5217 !> Soubroutine reads t_surf and t_wall data from restart file(s)
    5218 !kanani: Renamed this routine according to corresponging routines in PALM
    5219 !kanani: Modified the routine to match rrd_global, from where usm_rrd_local
    5220 !        shall be called in the future. This part has not been tested yet. (see virtual_flight_mod)
    5221 !        Also, I had some trouble with the allocation of t_surf, since this is a pointer.
    5222 !        So, I added some directives here.
     5338!> Soubroutine reads t_surf and t_wall data from restart files
    52235339!------------------------------------------------------------------------------!
    52245340    SUBROUTINE usm_rrd_local( i, k, nxlf, nxlc, nxl_on_file, nxrf, nxrc,       &
     
    66236739              surf_usm_h%albedo(:,m) = surface_params(ialbedo,ip)
    66246740           ENDIF
     6741!--        Albedo type is 0 (custom), others are replaced later
     6742           surf_usm_h%albedo_type(:,m) = 0
    66256743!--        Transmissivity
    66266744           IF ( surf_usm_h%transmissivity(m) < 0.0_wp )  THEN
     
    66366754           surf_usm_h%lambda_surf_green(m)  = surface_params(ilambdas,ip)
    66376755!           
    6638 !--        roughness relative to concrete
     6756!--        roughness length for momentum, heat and humidity
    66396757           surf_usm_h%z0(m) = surface_params(irough,ip)
    6640 !           
     6758           surf_usm_h%z0h(m) = surface_params(iroughh,ip)
     6759           surf_usm_h%z0q(m) = surface_params(iroughh,ip)
     6760!
    66416761!--        Surface skin layer heat capacity (J m−2 K−1 )
    66426762           surf_usm_h%c_surface(m) = surface_params(icsurf,ip)
     
    66766796              j  = surf_usm_v(l)%j(m)
    66776797              kw = surf_usm_v(l)%k(m)
    6678 
     6798             
    66796799              IF ( l == 3 )  THEN ! westward facing
    66806800                 iw = i
     
    66996819              ENDIF
    67006820
    6701               IF ( kw <= usm_par(ii,jw,iw) )  THEN
    6702 !--              pedestrant zone
     6821              IF ( iw < 0 .OR. jw < 0 ) THEN
     6822!--              wall on west or south border of the domain - assign default category
     6823                 IF ( kw <= roof_height_limit ) THEN
     6824                     surf_usm_v(l)%surface_types(m) = wall_category   !< default category for wall surface in wall zone
     6825                 ELSE
     6826                     surf_usm_v(l)%surface_types(m) = roof_category   !< default category for wall surface in roof zone
     6827                 END IF
     6828                 surf_usm_v(l)%albedo(:,m)    = -1.0_wp
     6829                 surf_usm_v(l)%thickness_wall(m) = -1.0_wp
     6830              ELSE IF ( kw <= usm_par(ii,jw,iw) )  THEN
     6831!--                 pedestrian zone
    67036832                 IF ( usm_par(ii+1,jw,iw) == 0 )  THEN
    6704                      surf_usm_v(l)%surface_types(m)  = pedestrant_category   !< default category for wall surface in pedestrant zone
     6833                     surf_usm_v(l)%surface_types(m)  = pedestrian_category   !< default category for wall surface in pedestrian zone
    67056834                     surf_usm_v(l)%albedo(:,m)    = -1.0_wp
    67066835                     surf_usm_v(l)%thickness_wall(m) = -1.0_wp
     
    67516880                 ENDIF
    67526881              ELSE
    6753 !--              something wrong
    6754                  CALL message( 'usm_read_urban_surface', 'PA0505', 1, 2, 0, 6, 0 )
     6882!
     6883!--              supply the default category
     6884                 IF ( kw <= roof_height_limit ) THEN
     6885                     surf_usm_v(l)%surface_types(m) = wall_category   !< default category for wall surface in wall zone
     6886                 ELSE
     6887                     surf_usm_v(l)%surface_types(m) = roof_category   !< default category for wall surface in roof zone
     6888                 END IF
     6889                 surf_usm_v(l)%albedo(:,m)    = -1.0_wp
     6890                 surf_usm_v(l)%thickness_wall(m) = -1.0_wp
    67556891              ENDIF
    6756 
    67576892!
    67586893!--           Find the type position
     
    67676902              IF ( ip == -99999 )  THEN
    67686903!--              wall category not found
    6769                  WRITE (message_string, "(A,I5,A,3I5)") 'wall category ', it,  &
     6904                 WRITE (message_string, "(A,I7,A,3I5)") 'wall category ', it,  &
    67706905                                        ' not found  for i,j,k=', iw,jw,kw
    6771                  CALL message( 'usm_read_urban_surface', 'PA0506', 1, 2, 0, 6, 0 )
     6906                 WRITE(9,*) message_string
    67726907              ENDIF
    67736908!
     
    67766911                 surf_usm_v(l)%albedo(:,m) = surface_params(ialbedo,ip)
    67776912              ENDIF
     6913!--           Albedo type is 0 (custom), others are replaced later
     6914              surf_usm_v(l)%albedo_type(:,m) = 0
    67786915!--           Transmissivity of the windows
    67796916              IF ( surf_usm_v(l)%transmissivity(m) < 0.0_wp )  THEN
     
    67846921              surf_usm_v(l)%emissivity(:,m) = surface_params(iemiss,ip)
    67856922!           
    6786 !--           heat conductivity λS between air and wall ( W m−2 K−1 )
     6923!--           heat conductivity lambda S between air and wall ( W m-2 K-1 )
    67876924              surf_usm_v(l)%lambda_surf(m) = surface_params(ilambdas,ip)
    67886925              surf_usm_v(l)%lambda_surf_window(m) = surface_params(ilambdas,ip)
    67896926              surf_usm_v(l)%lambda_surf_green(m) = surface_params(ilambdas,ip)
    67906927!           
    6791 !--           roughness relative to concrete
     6928!--           roughness length
    67926929              surf_usm_v(l)%z0(m) = surface_params(irough,ip)
     6930              surf_usm_v(l)%z0h(m) = surface_params(iroughh,ip)
     6931              surf_usm_v(l)%z0q(m) = surface_params(iroughh,ip)
    67936932!           
    6794 !--           Surface skin layer heat capacity (J m−2 K−1 )
     6933!--           Surface skin layer heat capacity (J m-2 K-1 )
    67956934              surf_usm_v(l)%c_surface(m) = surface_params(icsurf,ip)
    67966935              surf_usm_v(l)%c_surface_window(m) = surface_params(icsurf,ip)
     
    68096948                   surf_usm_v(l)%thickness_green(m) = surface_params(ithick,ip)
    68106949              ENDIF
    6811 !           
    6812 !--           volumetric heat capacity rho*C of the wall ( J m−3 K−1 )
     6950!
     6951!--           volumetric heat capacity rho*C of the wall ( J m-3 K-1 )
    68136952              surf_usm_v(l)%rho_c_wall(:,m) = surface_params(irhoC,ip)
    68146953              surf_usm_v(l)%rho_c_window(:,m) = surface_params(irhoC,ip)
    68156954              surf_usm_v(l)%rho_c_green(:,m) = surface_params(irhoC,ip)
    68166955!           
    6817 !--           thermal conductivity λH of the wall (W m−1 K−1 )
     6956!--           thermal conductivity lambda H of the wall (W m-1 K-1 )
    68186957              surf_usm_v(l)%lambda_h(:,m) = surface_params(ilambdah,ip)
    68196958              surf_usm_v(l)%lambda_h_window(:,m) = surface_params(ilambdah,ip)
     
    68486987   
    68496988    END SUBROUTINE usm_read_urban_surface_types
     6989
     6990
     6991!------------------------------------------------------------------------------!
     6992! Description:
     6993! ------------
     6994!
     6995!> This function advances through the list of local surfaces to find given
     6996!> x, y, d, z coordinates
     6997!------------------------------------------------------------------------------!
     6998    PURE FUNCTION advance_surface(isurfl_start, isurfl_stop, x, y, z, d) &
     6999            result(isurfl)
     7000
     7001        INTEGER(iwp), INTENT(in)                :: isurfl_start, isurfl_stop
     7002        INTEGER(iwp), INTENT(in)                :: x, y, z, d
     7003        INTEGER(iwp)                            :: isx, isy, isz, isd
     7004        INTEGER(iwp)                            :: isurfl
     7005
     7006        DO isurfl = isurfl_start, isurfl_stop
     7007            isx = surfl(ix, isurfl)
     7008            isy = surfl(iy, isurfl)
     7009            isz = surfl(iz, isurfl)
     7010            isd = surfl(id, isurfl)
     7011            IF ( isx==x .and. isy==y .and. isz==z .and. isd==d )  RETURN
     7012        ENDDO
     7013
     7014!--     coordinate not found
     7015        isurfl = -1
     7016
     7017    END FUNCTION
     7018
     7019
     7020!------------------------------------------------------------------------------!
     7021! Description:
     7022! ------------
     7023!
     7024!> This subroutine reads temperatures of respective material layers in walls,
     7025!> roofs and ground from input files. Data in the input file must be in
     7026!> standard order, i.e. horizontal surfaces first ordered by x, y and then
     7027!> vertical surfaces ordered by x, y, direction, z
     7028!------------------------------------------------------------------------------!
     7029    SUBROUTINE usm_read_wall_temperature
     7030
     7031        INTEGER(iwp)                                          :: i, j, k, d, ii, iline
     7032        INTEGER(iwp)                                          :: isurfl
     7033        REAL(wp)                                              :: rtsurf
     7034        REAL(wp), DIMENSION(nzb_wall:nzt_wall+1)              :: rtwall
     7035
     7036
     7037
     7038
     7039        DO  ii = 0, io_blocks-1
     7040            IF ( ii == io_group )  THEN
     7041
     7042!--             open wall temperature file
     7043                OPEN( 152, file='WALL_TEMPERATURE'//coupling_char, action='read', &
     7044                           status='old', form='formatted', err=15 )
     7045
     7046                isurfl = 0
     7047                iline = 1
     7048                DO
     7049                    rtwall = -9999.0_wp  !< for incomplete lines
     7050                    READ( 152, *, err=13, end=14 )  i, j, k, d, rtsurf, rtwall
     7051
     7052                    IF ( nxl <= i .and. i <= nxr .and. &
     7053                        nys <= j .and. j <= nyn)  THEN  !< local processor
     7054!--                     identify surface id
     7055                        isurfl = advance_surface(isurfl+1, nsurfl, i, j, k, d)
     7056                        IF ( isurfl == -1 )  THEN
     7057                            WRITE(message_string, '(a,4i5,a,i5,a)') 'Coordinates (xyzd) ', i, j, k, d, &
     7058                                ' on line ', iline, &
     7059                                ' in file WALL_TEMPERATURE are either not present or out of standard order of surfaces.'
     7060                            CALL message( 'usm_read_wall_temperature', 'PA0521', 1, 2, 0, 6, 0 )
     7061                        ENDIF
     7062
     7063!--                     assign temperatures
     7064                        IF ( d == 0 ) THEN
     7065                           t_surf_h(isurfl) = rtsurf
     7066                           t_wall_h(:,isurfl) = rtwall(:)
     7067                        ELSE
     7068                           t_surf_v(d-1)%t(isurfl) = rtsurf
     7069                           t_wall_v(d-1)%t(:,isurfl) = rtwall(:)
     7070                        ENDIF
     7071                    ENDIF
     7072
     7073                    iline = iline + 1
     7074                    CYCLE
     7075 13                 WRITE(message_string, '(a,i5,a)') 'Error reading line ', iline, &
     7076                        ' in file WALL_TEMPERATURE.'
     7077                    CALL message( 'usm_read_wall_temperature', 'PA0522', 1, 2, 0, 6, 0 )
     7078                ENDDO
     7079 14             CLOSE(152)
     7080                CYCLE
     7081 15             message_string = 'file WALL_TEMPERATURE'//TRIM(coupling_char)//' does not exist'
     7082                CALL message( 'usm_read_wall_temperature', 'PA0523', 1, 2, 0, 6, 0 )
     7083            ENDIF
     7084#if defined( __parallel ) && ! defined ( __check )
     7085            CALL MPI_BARRIER( comm2d, ierr )
     7086#endif
     7087        ENDDO
     7088
     7089        CALL location_message( '    wall layer temperatures read', .TRUE. )
     7090
     7091    END SUBROUTINE usm_read_wall_temperature
     7092
    68507093
    68517094
     
    69457188           IF ( surf_usm_h%r_a_window(m) < 1.0_wp )                            &
    69467189               surf_usm_h%r_a_window(m) = 1.0_wp
    6947 
    6948                
    6949 !--        the parameterization is developed originally for larger scales
    6950 !--        (compare with remark in TUF-3D)
    6951 !--        our first experiences show that the parameterization underestimates
    6952 !--        r_a in meter resolution.
    6953 !--        A temporary solution would be multiplication by magic constant :-(.
    6954 !--        For the moment this is comment out.
    6955            surf_usm_h%r_a(m)        = surf_usm_h%r_a(m)        !* ra_horiz_coef
    6956            surf_usm_h%r_a_window(m) = surf_usm_h%r_a_window(m) !* ra_horiz_coef
    6957            surf_usm_h%r_a_green(m)  = surf_usm_h%r_a_green(m)  !* ra_horiz_coef
    69587190               
    69597191!--        factor for shf_eb
     
    70707302!--        calculate fluxes
    70717303!--        rad_net_l is never used!           
    7072            surf_usm_h%rad_net_l(m) = surf_usm_h%frac(0,m) *                         &
    7073                                      ( surf_usm_h%rad_net_l(m) +                    &
    7074                                      3.0_wp * sigma_sb *                            &
    7075                                      t_surf_h(m)**4 - 4.0_wp * sigma_sb *           &
    7076                                      t_surf_h(m)**3 * t_surf_h_p(m) )               &
    7077                                    + surf_usm_h%frac(2,m) *                         &
    7078                                      ( surf_usm_h%rad_net_l(m) +                    &
    7079                                      3.0_wp * sigma_sb *                            &
    7080                                      t_surf_window_h(m)**4 - 4.0_wp * sigma_sb *    &
    7081                                      t_surf_window_h(m)**3 * t_surf_window_h_p(m) ) &
    7082                                    + surf_usm_h%frac(1,m) *                         &
    7083                                      ( surf_usm_h%rad_net_l(m) +                    &
    7084                                      3.0_wp * sigma_sb *                            &
    7085                                      t_surf_green_h(m)**4 - 4.0_wp * sigma_sb *     &
    7086                                      t_surf_green_h(m)**3 * t_surf_green_h_p(m) )
     7304           surf_usm_h%rad_net_l(m) = surf_usm_h%rad_net_l(m) +                           &
     7305                                     surf_usm_h%frac(0,m) *                              &
     7306                                     sigma_sb * surf_usm_h%emissivity(0,m) *             &
     7307                                     ( t_surf_h_p(m)**4 - t_surf_h(m)**4 )               &
     7308                                    + surf_usm_h%frac(2,m) *                         &
     7309                                     sigma_sb * surf_usm_h%emissivity(2,m) *             &
     7310                                     ( t_surf_window_h_p(m)**4 - t_surf_window_h(m)**4 ) &
     7311                                    + surf_usm_h%frac(1,m) *                         &
     7312                                     sigma_sb * surf_usm_h%emissivity(1,m) *             &
     7313                                     ( t_surf_green_h_p(m)**4 - t_surf_green_h(m)**4 )
     7314
    70877315           surf_usm_h%wghf_eb(m)   = lambda_surface *                         &
    70887316                                      ( t_surf_h_p(m) - t_wall_h(nzb_wall,m) )
     
    72697497
    72707498!--           calculate fluxes
    7271 !--           rad_net_l is never used!           
     7499!--           prognostic rad_net_l is used just for output!           
    72727500              surf_usm_v(l)%rad_net_l(m) = surf_usm_v(l)%frac(0,m) *                             &
    72737501                                           ( surf_usm_v(l)%rad_net_l(m) +                            &
     
    73227550           dtime = mod(simulated_time + time_utc_init, 24.0_wp*3600.0_wp)
    73237551           dhour = INT(dtime/3600.0_wp)
    7324 !--        linear interpolation of coeficient
    7325            acoef = (REAL(dhour+1,wp)-dtime/3600.0_wp)*aheatprof(dhour) + (dtime/3600.0_wp-REAL(dhour,wp))*aheatprof(dhour+1)
    7326 
    7327            DO m = 1, surf_usm_h%ns
    7328 !
     7552           DO m = 1, naheatlayers
    73297553!--           Get indices of respective grid point
    73307554              i = surf_usm_h%i(m)
    73317555              j = surf_usm_h%j(m)
    73327556              k = surf_usm_h%k(m)
    7333 
    7334               IF ( aheat(j,i) > 0.0_wp )  THEN
    7335 !--              TODO the increase of pt in box i,j,nzb_s_inner(j,i)+1 in time dt_3d
     7557              IF ( k > get_topography_top_index_ji( j, i, 's' )  .AND. &
     7558                   k <= naheatlayers )  THEN
     7559!--              increase of pt in box i,j,k in time dt_3d
    73367560!--              given to anthropogenic heat aheat*acoef (W*m-2)
    7337 !--              k = nzb_s_inner(j,i)+1
    7338 !--              pt(k,j,i) = pt(k,j,i) + aheat(j,i)*acoef*dt_3d/(exn(k)*rho_cp*dz)
    7339 !--              Instead of this, we can adjust shf in case AH only at surface
    7340                  surf_usm_h%shf(m) = surf_usm_h%shf(m) +                       &
    7341                                    aheat(j,i) * acoef * ddx * ddy / cp
     7561!--              linear interpolation of coeficient
     7562                 acoef = (REAL(dhour+1,wp)-dtime/3600.0_wp)*aheatprof(k,dhour) + &
     7563                         (dtime/3600.0_wp-REAL(dhour,wp))*aheatprof(k,dhour+1)
     7564                 IF ( aheat(k,j,i) > 0.0_wp )  THEN
     7565                    pt(k,j,i) = pt(k,j,i) + aheat(k,j,i)*acoef*dt_3d/(exn(k)*rho_cp*dz)
     7566                 ENDIF
    73427567              ENDIF
    73437568           ENDDO
     
    74377662
    74387663!------------------------------------------------------------------------------!
    7439 !
    74407664! Description:
    74417665! ------------
    7442 !> Subroutine writes t_surf and t_wall data into restart files. It is necessary
    7443 !> to write start_index and end_index several times
    7444 !kanani: Renamed this routine according to corresponging routines in PALM
    7445 !kanani: Modified the routine to match wrd_global, from where usm_wrd_local
    7446 !        shall be called in the future. This part has not been tested yet. (see virtual_flight_mod)
    7447 !        Also, I had some trouble with the allocation of t_surf, since this is a pointer.
    7448 !        So, I added some directives here.
     7666!> Subroutine writes t_surf and t_wall data into restart files
    74497667!------------------------------------------------------------------------------!
    74507668    SUBROUTINE usm_wrd_local
Note: See TracChangeset for help on using the changeset viewer.