Ignore:
Timestamp:
May 30, 2017 5:47:52 PM (4 years ago)
Author:
suehring
Message:

Adjustments according new topography and surface-modelling concept implemented

File:
1 edited

Legend:

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

    r2214 r2232  
    2121! Current revisions:
    2222! ------------------
    23 !
     23! Adjustments according to new surface-type structure. Remove usm_wall_heat_flux;
     24! insteat, heat fluxes are directly applied in diffusion_s.
    2425!
    2526! Former revisions:
     
    103104!>    reflection step using MPI_Alltoall instead of current MPI_Allgather.
    104105!>
     106!> 3. Temporarily large values of surface heat flux can be observed, up to
     107!>    1.2 Km/s, which seem to be not realistic.
     108!>
    105109!> @todo Check optimizations for RMA operations
    106110!> @todo Alternatives for MPI_WIN_ALLOCATE? (causes problems with openmpi)
     
    111115
    112116    USE arrays_3d,                                                             &
    113         ONLY:  zu, pt, pt_1, pt_2, p, ol, shf, ts, us, u, v, w, hyp, tend
     117        ONLY:  zu, pt, pt_1, pt_2, p, u, v, w, hyp, tend
    114118
    115119    USE cloud_parameters,                                                      &
     
    137141    USE indices,                                                               &
    138142        ONLY:  nx, ny, nnx, nny, nnz, nxl, nxlg, nxr, nxrg, nyn, nyng, nys,    &
    139                nysg, nzb_s_inner, nzb_s_outer, nzb, nzt, nbgp
     143               nysg, nzb, nzt, nbgp, wall_flags_0
    140144
    141145    USE, INTRINSIC :: iso_c_binding
     
    157161    USE statistics,                                                            &
    158162        ONLY:  hom, statistic_regions
     163
     164    USE surface_mod
    159165
    160166               
     
    265271        REAL(wp)                                   :: rtransp           !<
    266272    END TYPE
     273!
     274!-- Type for surface temperatures at vertical walls. Is not necessary for horizontal walls.
     275    TYPE t_surf_vertical
     276       REAL(wp), DIMENSION(:), ALLOCATABLE         :: t
     277    END TYPE t_surf_vertical
     278!
     279!-- Type for wall temperatures at vertical walls. Is not necessary for horizontal walls.
     280    TYPE t_wall_vertical
     281       REAL(wp), DIMENSION(:,:), ALLOCATABLE       :: t
     282    END TYPE t_wall_vertical
    267283
    268284!-- arrays for calculation of svf and csf
     
    320336    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfins_av       !< average of array of residua of sw radiation absorbed in surface after last reflection
    321337    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinl_av       !< average of array of residua of lw radiation absorbed in surface after last reflection
    322     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfhf_av        !< average of total radiation flux incoming to minus outgoing from local surface   
     338    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfhf_av        !< average of total radiation flux incoming to minus outgoing from local surface 
    323339   
    324340!-- block variables needed for calculation of the plant canopy model inside the urban surface model
     
    367383!-- surface and material model variables for walls, ground, roofs
    368384!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    369     INTEGER(iwp), DIMENSION(:), ALLOCATABLE        :: surface_types      !< array of types of wall parameters
    370385    REAL(wp), DIMENSION(:), ALLOCATABLE            :: zwn                !< normalized wall layer depths (m)
    371     REAL(wp), DIMENSION(:,:), ALLOCATABLE          :: ddz_wall           !< 1/dz_wall
    372     REAL(wp), DIMENSION(:,:), ALLOCATABLE          :: ddz_wall_stag      !< 1/dz_wall_stag
    373     REAL(wp), DIMENSION(:,:), ALLOCATABLE          :: dz_wall            !< wall grid spacing (center-center)
    374     REAL(wp), DIMENSION(:,:), ALLOCATABLE          :: dz_wall_stag       !< wall grid spacing (edge-edge)
    375     REAL(wp), DIMENSION(:,:), ALLOCATABLE          :: zw                 !< wall layer depths (m)
    376386
    377387#if defined( __nopointer )
    378     REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf             !< wall surface temperature (K)
    379     REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_p           !< progn. wall surface temperature (K)
     388    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_h           !< wall surface temperature (K) at horizontal walls
     389    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_h_p         !< progn. wall surface temperature (K) at horizontal walls
     390
     391    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  ::  t_surf_v
     392    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  ::  t_surf_v_p
    380393#else
    381     REAL(wp), DIMENSION(:), POINTER                :: t_surf
    382     REAL(wp), DIMENSION(:), POINTER                :: t_surf_p
    383 
    384     REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_1
    385     REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_2
     394    REAL(wp), DIMENSION(:), POINTER                :: t_surf_h
     395    REAL(wp), DIMENSION(:), POINTER                :: t_surf_h_p
     396
     397    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_h_1
     398    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_h_2
     399
     400    TYPE(t_surf_vertical), DIMENSION(:), POINTER ::  t_surf_v
     401    TYPE(t_surf_vertical), DIMENSION(:), POINTER ::  t_surf_v_p
     402
     403    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  :: t_surf_v_1
     404    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  :: t_surf_v_2
    386405#endif
    387406    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_av          !< average of wall surface temperature (K)
     
    394413!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    395414!-- parameters of the land, roof and wall surfaces
    396     LOGICAL,  DIMENSION(:), ALLOCATABLE            :: isroof_surf        !< is the surface the part of a roof
    397415    REAL(wp), DIMENSION(:), ALLOCATABLE            :: albedo_surf        !< albedo of the surface
    398416!-- parameters of the wall surfaces
    399     REAL(wp), DIMENSION(:), ALLOCATABLE            :: c_surface          !< heat capacity of the wall surface skin ( J m−2 K−1 )
    400417    REAL(wp), DIMENSION(:), ALLOCATABLE            :: emiss_surf         !< emissivity of the wall surface
    401     REAL(wp), DIMENSION(:), ALLOCATABLE            :: lambda_surf        !< heat conductivity λS between air and surface ( W m−2 K−1 )
    402    
    403 !-- parameters of the walls material
    404     REAL(wp), DIMENSION(:), ALLOCATABLE            :: thickness_wall     !< thickness of the wall, roof and soil layers
    405     REAL(wp), DIMENSION(:,:), ALLOCATABLE          :: rho_c_wall         !< volumetric heat capacity of the material ( J m-3 K-1 ) (= 2.19E6)
    406     REAL(wp), DIMENSION(:,:), ALLOCATABLE          :: lambda_h           !< heat conductivity λT of the material ( W m-1 K-1 )
    407     REAL(wp), DIMENSION(:), ALLOCATABLE            :: roughness_wall     !< roughness relative to concrete
    408    
    409 !-- output wall heat flux arrays
    410     REAL(wp), DIMENSION(:), ALLOCATABLE            :: wshf               !< kinematic wall heat flux of sensible heat (needed for diffusion_s!<)
    411     REAL(wp), DIMENSION(:), ALLOCATABLE            :: wshf_eb            !< wall heat flux of sensible heat in wall normal direction
    412     REAL(wp), DIMENSION(:), ALLOCATABLE            :: wshf_eb_av         !< average of wshf_eb
    413     REAL(wp), DIMENSION(:), ALLOCATABLE            :: wghf_eb            !< wall ground heat flux
    414     REAL(wp), DIMENSION(:), ALLOCATABLE            :: wghf_eb_av         !< average of wghf_eb
    415418
    416419#if defined( __nopointer )
    417     REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET  :: t_wall             !< Wall temperature (K)
    418     REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET  :: t_wall_av          !< Average of t_wall
    419     REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET  :: t_wall_p           !< Prog. wall temperature (K)
     420    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_wall_h             !< Wall temperature (K)
     421    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_wall_h_av          !< Average of t_wall
     422    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_wall_h_p           !< Prog. wall temperature (K)
     423
     424    TYPE(t_wall_vertical), DIMENSION(0:3), TARGET  :: t_wall_v             !< Wall temperature (K)
     425    TYPE(t_wall_vertical), DIMENSION(0:3), TARGET  :: t_wall_v_av          !< Average of t_wall
     426    TYPE(t_wall_vertical), DIMENSION(0:3), TARGET  :: t_wall_v_p           !< Prog. wall temperature (K)
    420427#else
    421     REAL(wp), DIMENSION(:,:), POINTER              :: t_wall, t_wall_p
    422     REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET  :: t_wall_av, t_wall_1, t_wall_2
     428    REAL(wp), DIMENSION(:,:), POINTER                :: t_wall_h, t_wall_h_p
     429    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET    :: t_wall_h_av, t_wall_h_1, t_wall_h_2
     430
     431    TYPE(t_wall_vertical), DIMENSION(:), POINTER   :: t_wall_v, t_wall_v_p
     432    TYPE(t_wall_vertical), DIMENSION(0:3), TARGET  :: t_wall_v_av, t_wall_v_1, t_wall_v_2
    423433#endif
    424434
     
    488498       MODULE PROCEDURE usm_swap_timelevel
    489499    END INTERFACE usm_swap_timelevel
    490    
    491     INTERFACE usm_wall_heat_flux
    492        MODULE PROCEDURE usm_wall_heat_flux
    493        MODULE PROCEDURE usm_wall_heat_flux_ij
    494     END INTERFACE usm_wall_heat_flux
    495    
     500       
    496501    INTERFACE usm_write_restart_data
    497502       MODULE PROCEDURE usm_write_restart_data
     
    508513           usm_energy_balance_land, usm_energy_balance_wall, nrefsteps,        &
    509514           usm_init_urban_surface, usm_radiation, usm_read_restart_data,       &
    510            usm_wall_heat_flux,                                                 &
    511515           usm_surface_energy_balance, usm_material_heat_model,                &
    512516           usm_swap_timelevel, usm_check_data_output, usm_average_3d_data,     &
     
    531535        IMPLICIT NONE
    532536       
    533         INTEGER(iwp)                            :: i, j, k, d, l, ir, jr, ids
    534         INTEGER(iwp)                            :: nzubl, nzutl, isurf, ipcgb
    535         INTEGER(iwp)                            :: procid
     537        INTEGER(iwp) :: i, j, k, d, l, ir, jr, ids, m
     538        INTEGER(iwp) :: k_topo     !< vertical index indicating topography top for given (j,i)
     539        INTEGER(iwp) :: k_topo2    !< vertical index indicating topography top for given (j,i)
     540        INTEGER(iwp) :: nzubl, nzutl, isurf, ipcgb
     541        INTEGER(iwp) :: procid
    536542
    537543       
     
    543549        CALL location_message( '', .TRUE. )
    544550        CALL location_message( '    allocation of needed arrays', .TRUE. )
    545 !--     find nzub, nzut, nzu
    546         nzubl = minval(nzb_s_inner(nys:nyn,nxl:nxr))
    547         nzutl = maxval(nzb_s_inner(nys:nyn,nxl:nxr))
     551!
     552!--     Find nzub, nzut, nzu via wall_flag_0 array (nzb_s_inner will be
     553!--     removed later). The following contruct finds the lowest / largest index
     554!--     for any upward-facing wall (see bit 12).
     555        nzubl = MINVAL(                                                        &
     556                    MAXLOC(                                                    &
     557                          MERGE( 1, 0,                                         &
     558                                 BTEST( wall_flags_0(:,nys:nyn,nxl:nxr), 12 )  &
     559                               ), DIM = 1                                      &
     560                          ) - 1                                                &
     561                            )
     562        nzutl = MAXVAL(                                                        &
     563                   MAXLOC(                                                     &
     564                          MERGE( 1, 0,                                         &
     565                                 BTEST( wall_flags_0(:,nys:nyn,nxl:nxr), 12 )  &
     566                               ), DIM = 1                                      &
     567                          ) - 1                                                &
     568                            )
    548569        nzubl = max(nzubl,nzb)
     570
    549571       
    550572        IF ( plant_canopy )  THEN
     
    559581            DO i = nxl, nxr
    560582                DO j = nys, nyn
     583!
     584!--                 Find topography top index
     585                    k_topo = MAXLOC( MERGE( 1, 0,                              &
     586                                             BTEST( wall_flags_0(:,j,i), 12 )  &
     587                                           ), DIM = 1                          &
     588                                   ) - 1
    561589                    DO k = nzt+1, 0, -1
    562590                        IF ( lad_s(k,j,i) /= 0.0_wp )  THEN
    563591!--                         we are at the top of the pcs
    564                             pct(j,i) = k + nzb_s_inner(j,i)
     592                            pct(j,i) = k + k_topo
    565593                            pch(j,i) = k
    566594                            npcbl = npcbl + pch(j,i)
     
    601629        CALL location_message( '    calculation of indices for surfaces', .TRUE. )
    602630        nsurfl = 0
    603 !--     calculate land surface and roof
    604         startland = nsurfl+1
    605         nsurfl = nsurfl+(nxr-nxl+1)*(nyn-nys+1)
    606         endland = nsurfl
    607         nlands = endland-startland+1
    608 
    609 !--     calculation of the walls
     631!
     632!--     Number of land- and roof surfaces. Note, since horizontal surface elements
     633!--     are already counted in surface_mod, in case be simply reused here.
     634        startland = 1
     635        nsurfl    = surf_usm_h%ns
     636        endland   = nsurfl
     637        nlands    = endland-startland+1
     638
     639!
     640!--     Number of vertical surfaces. As vertical surfaces are already
     641!--     counted in surface mod, it can be reused here.
    610642        startwall = nsurfl+1
    611         DO i = nxl, nxr
    612             DO j = nys, nyn
    613 !--             test for walls
    614 !--             (we don't use array flags because it isn't calculated in case of masking_method=.T.)
    615                 DO ids = 1, 4  !-- four wall directions
    616                     jr = min(max(j-jdir(ids),0),ny)
    617                     ir = min(max(i-idir(ids),0),nx)
    618                     nsurfl = nsurfl + max(0, nzb_s_inner(jr,ir)-nzb_s_inner(j,i))
    619                 ENDDO
    620             ENDDO
    621         ENDDO
     643        nsurfl = nsurfl + surf_usm_v(0)%ns + surf_usm_v(1)%ns +        &
     644                          surf_usm_v(2)%ns + surf_usm_v(3)%ns
    622645        endwall = nsurfl
    623646        nwalls = endwall-startwall+1
    624        
    625 !--     range of energy balance surfaces
     647
     648       
     649!--     range of energy balance surfaces  ! will be treated separately by surf_usm_h and surf_usm_v
    626650        nenergy = 0
    627651        IF ( usm_energy_balance_land )  THEN
     
    641665!--     block of virtual surfaces
    642666!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    643 !--     calculate sky surfaces
     667!--     calculate sky surfaces  ! not used so far!
    644668        startsky = nsurfl+1
    645669        nsurfl = nsurfl+(nxr-nxl+1)*(nyn-nys+1)
     
    656680        ijdb = RESHAPE( (/ nxl,nxr,nyn,nyn,nxl,nxr,nys,nys,nxr,nxr,nys,nyn,nxl,nxl,nys,nyn /), (/4, 4/) )
    657681!--     calulation of the free borders of the domain
    658         DO ids = 6,9
    659             IF ( isborder(ids) )  THEN
    660 !--             free border of the domain in direction ids
    661                 DO i = ijdb(1,ids), ijdb(2,ids)
    662                     DO j = ijdb(3,ids), ijdb(4,ids)
    663                         k = nzut - max(nzb_s_inner(j,i), nzb_s_inner(j-jdir(ids),i-idir(ids)))
    664                         nsurfl = nsurfl + k
    665                     ENDDO
    666                 ENDDO
    667             ENDIF
     682        DO  ids = 6,9
     683           IF ( isborder(ids) )  THEN
     684!--           free border of the domain in direction ids
     685              DO  i = ijdb(1,ids), ijdb(2,ids)
     686                 DO  j = ijdb(3,ids), ijdb(4,ids)
     687
     688                    k_topo =   MAXLOC( MERGE( 1, 0,                            &
     689                                             BTEST( wall_flags_0(:,j,i), 12 )  &
     690                                            ), DIM = 1                         &
     691                                     ) - 1
     692                    k_topo2 =  MAXLOC( MERGE( 1, 0,                            &
     693                                             BTEST( wall_flags_0(:,j-jdir(ids),i-idir(ids)), 12 )  &
     694                                            ), DIM = 1                         &
     695                                     ) - 1
     696
     697                    k = nzut - MAX( k_topo, k_topo2 )
     698                    nsurfl = nsurfl + k
     699                 ENDDO
     700              ENDDO
     701           ENDIF
    668702        ENDDO
    669703       
     
    676710            DO i = nxl, nxr
    677711                DO j = nys, nyn
    678                     DO k = nzb_s_inner(j,i)+1, pct(j,i)
     712!
     713!--                 Find topography top index
     714                    k_topo = MAXLOC( MERGE( 1, 0,                              &
     715                                             BTEST( wall_flags_0(:,j,i), 12 )  &
     716                                          ), DIM = 1                           &
     717                                   ) - 1
     718                    DO k = k_topo + 1, pct(j,i)
    679719                        ipcgb = ipcgb + 1
    680720                        gridpcbl(k,j,i) = ipcgb
     
    689729
    690730!--     fill surfl
    691         ALLOCATE(surfl(4,nsurfl))
     731        ALLOCATE(surfl(5,nsurfl))
    692732        isurf = 0
    693733       
     
    695735        DO i = nxl, nxr
    696736            DO j = nys, nyn
    697                 isurf = isurf + 1
    698                 k = nzb_s_inner(j,i)+1
    699                 surfl(:,isurf) = (/iroof,k,j,i/)
     737               DO  m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
     738                  k = surf_usm_h%k(m)
     739
     740                  isurf = isurf + 1
     741                  surfl(:,isurf) = (/iroof,k,j,i,m/)
     742               ENDDO
    700743            ENDDO
    701744        ENDDO
     
    704747        DO i = nxl, nxr
    705748            DO j = nys, nyn
    706                 DO ids = 1, 4  !> four wall directions
    707                     jr = min(max(j-jdir(ids),0),ny)
    708                     ir = min(max(i-idir(ids),0),nx)
    709                     DO k = nzb_s_inner(j,i)+1, nzb_s_inner(jr,ir)
    710                         isurf = isurf + 1
    711                         surfl(:,isurf) = (/ids,k,j,i/)
    712                     ENDDO
    713                 ENDDO
     749               l = 0
     750               DO  m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i)
     751                  k = surf_usm_v(l)%k(m)
     752
     753                  isurf          = isurf + 1
     754                  surfl(:,isurf) = (/2,k,j,i,m/)
     755               ENDDO
     756               l = 1
     757               DO  m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i)
     758                  k = surf_usm_v(l)%k(m)
     759
     760                  isurf          = isurf + 1
     761                  surfl(:,isurf) = (/1,k,j,i,m/)
     762               ENDDO
     763               l = 2
     764               DO  m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i)
     765                  k = surf_usm_v(l)%k(m)
     766
     767                  isurf          = isurf + 1
     768                  surfl(:,isurf) = (/4,k,j,i,m/)
     769               ENDDO
     770               l = 3
     771               DO  m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i)
     772                  k = surf_usm_v(l)%k(m)
     773
     774                  isurf          = isurf + 1
     775                  surfl(:,isurf) = (/3,k,j,i,m/)
     776               ENDDO
    714777            ENDDO
    715778        ENDDO
     
    720783                isurf = isurf + 1
    721784                k = nzut
    722                 surfl(:,isurf) = (/isky,k,j,i/)
     785                surfl(:,isurf) = (/isky,k,j,i,-1/)
    723786            ENDDO
    724787        ENDDO
     
    730793                DO i = ijdb(1,ids), ijdb(2,ids)
    731794                    DO j = ijdb(3,ids), ijdb(4,ids)
    732                         DO k = max(nzb_s_inner(j,i),nzb_s_inner(j-jdir(ids),i-idir(ids)))+1, nzut
     795                        k_topo =   MAXLOC( MERGE( 1, 0,                        &
     796                                             BTEST( wall_flags_0(:,j,i), 12 )  &
     797                                                ), DIM = 1                     &
     798                                         ) - 1
     799                        k_topo2 =  MAXLOC( MERGE( 1, 0,                        &
     800                                             BTEST( wall_flags_0(:,j-jdir(ids),i-idir(ids)), 12 )  &
     801                                                ), DIM = 1                     &
     802                                         ) - 1
     803                        DO k = MAX(k_topo,k_topo2)+1, nzut
    733804                            isurf = isurf + 1
    734                             surfl(:,isurf) = (/ids,k,j,i/)
     805                            surfl(:,isurf) = (/ids,k,j,i,-1/)
    735806                        ENDDO
    736807                    ENDDO
     
    755826        surfstart(numprocs) = k
    756827        nsurf = k
    757         ALLOCATE(surf(4,nsurf))
     828        ALLOCATE(surf(5,nsurf))
    758829       
    759830#if defined( __parallel )
    760         CALL MPI_AllGatherv(surfl, nsurfl*4, MPI_INTEGER, surf, nsurfs*4, surfstart*4, MPI_INTEGER, comm2d, ierr)
     831        CALL MPI_AllGatherv(surfl, nsurfl*5, MPI_INTEGER, surf, nsurfs*5, surfstart*5, MPI_INTEGER, comm2d, ierr)
    761832#else
    762833        surf = surfl
     
    787858        ALLOCATE( surfouts(nsurf) ) !TODO: global surfaces without virtual
    788859        ALLOCATE( surfoutl(nsurf) ) !TODO: global surfaces without virtual
    789         ALLOCATE( surfhf(startenergy:endenergy) )
    790         ALLOCATE( rad_net_l(startenergy:endenergy) )
     860
     861
     862
     863!
     864!--     Allocate radiation arrays which are part of the new data type.
     865!--     For horizontal surfaces.
     866        ALLOCATE( surf_usm_h%surfhf(1:surf_usm_h%ns)    )
     867        ALLOCATE( surf_usm_h%rad_net_l(1:surf_usm_h%ns) )
     868!
     869!--  New
     870        ALLOCATE( surf_usm_h%rad_in_sw(1:surf_usm_h%ns)  )
     871        ALLOCATE( surf_usm_h%rad_out_sw(1:surf_usm_h%ns) )
     872        ALLOCATE( surf_usm_h%rad_in_lw(1:surf_usm_h%ns)  )
     873        ALLOCATE( surf_usm_h%rad_out_lw(1:surf_usm_h%ns) )
     874!
     875!--     For vertical surfaces
     876        DO  l = 0, 3
     877           ALLOCATE( surf_usm_v(l)%surfhf(1:surf_usm_v(l)%ns)    )
     878           ALLOCATE( surf_usm_v(l)%rad_net_l(1:surf_usm_v(l)%ns) )
     879           ALLOCATE( surf_usm_v(l)%rad_in_sw(1:surf_usm_v(l)%ns)  )
     880           ALLOCATE( surf_usm_v(l)%rad_out_sw(1:surf_usm_v(l)%ns) )
     881           ALLOCATE( surf_usm_v(l)%rad_in_lw(1:surf_usm_v(l)%ns)  )
     882           ALLOCATE( surf_usm_v(l)%rad_out_lw(1:surf_usm_v(l)%ns) )
     883        ENDDO
    791884
    792885!--     Wall surface model
     
    794887       
    795888!--     allocate array of wall types and wall parameters
    796         ALLOCATE ( surface_types(startenergy:endenergy) )
     889        ALLOCATE ( surf_usm_h%surface_types(1:surf_usm_h%ns) )
     890        DO  l = 0, 3
     891           ALLOCATE( surf_usm_v(l)%surface_types(1:surf_usm_v(l)%ns) )
     892        ENDDO
    797893       
    798894!--     broadband albedo of the land, roof and wall surface
     
    802898        ALLOCATE ( albedo_surf(nsurfl) )
    803899        albedo_surf = 1.0_wp
     900        ALLOCATE ( surf_usm_h%albedo_surf(1:surf_usm_h%ns) )
     901        DO  l = 0, 3
     902           ALLOCATE( surf_usm_v(l)%albedo_surf(1:surf_usm_v(l)%ns) )
     903        ENDDO
    804904       
    805 !--     wall and roof surface parameters
    806         ALLOCATE ( isroof_surf(startenergy:endenergy) )
     905!--     wall and roof surface parameters. First for horizontal surfaces
    807906        ALLOCATE ( emiss_surf(startenergy:endenergy) )
    808         ALLOCATE ( lambda_surf(startenergy:endenergy) )
    809         ALLOCATE ( c_surface(startenergy:endenergy) )
    810         ALLOCATE ( roughness_wall(startenergy:endenergy) )
     907
     908        ALLOCATE ( surf_usm_h%isroof_surf(1:surf_usm_h%ns)    )
     909        ALLOCATE ( surf_usm_h%emiss_surf(1:surf_usm_h%ns)     )
     910        ALLOCATE ( surf_usm_h%lambda_surf(1:surf_usm_h%ns)    )
     911        ALLOCATE ( surf_usm_h%c_surface(1:surf_usm_h%ns)      )
     912        ALLOCATE ( surf_usm_h%roughness_wall(1:surf_usm_h%ns) )
     913!
     914!--     For vertical surfaces.
     915        DO  l = 0, 3
     916           ALLOCATE ( surf_usm_v(l)%emiss_surf(1:surf_usm_v(l)%ns)     )
     917           ALLOCATE ( surf_usm_v(l)%lambda_surf(1:surf_usm_v(l)%ns)    )
     918           ALLOCATE ( surf_usm_v(l)%c_surface(1:surf_usm_v(l)%ns)      )
     919           ALLOCATE ( surf_usm_v(l)%roughness_wall(1:surf_usm_v(l)%ns) )
     920        ENDDO
    811921       
    812 !--     allocate wall and roof material parameters
    813         ALLOCATE ( thickness_wall(startenergy:endenergy) )
    814         ALLOCATE ( lambda_h(nzb_wall:nzt_wall,startenergy:endenergy) )
    815         ALLOCATE ( rho_c_wall(nzb_wall:nzt_wall,startenergy:endenergy) )
    816 
    817 !--     allocate wall and roof layers sizes
     922!--     allocate wall and roof material parameters. First for horizontal surfaces
     923        ALLOCATE ( surf_usm_h%thickness_wall(1:surf_usm_h%ns)               )
     924        ALLOCATE ( surf_usm_h%lambda_h(nzb_wall:nzt_wall,1:surf_usm_h%ns)   )
     925        ALLOCATE ( surf_usm_h%rho_c_wall(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
     926!
     927!--     For vertical surfaces.
     928        DO  l = 0, 3
     929           ALLOCATE ( surf_usm_v(l)%thickness_wall(1:surf_usm_v(l)%ns)               )
     930           ALLOCATE ( surf_usm_v(l)%lambda_h(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns)   )
     931           ALLOCATE ( surf_usm_v(l)%rho_c_wall(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
     932        ENDDO
     933
     934!--     allocate wall and roof layers sizes. For horizontal surfaces.
    818935        ALLOCATE ( zwn(nzb_wall:nzt_wall) )
    819         ALLOCATE ( dz_wall(nzb_wall:nzt_wall+1, startenergy:endenergy) )
    820         ALLOCATE ( ddz_wall(nzb_wall:nzt_wall+1, startenergy:endenergy) )
    821         ALLOCATE ( dz_wall_stag(nzb_wall:nzt_wall, startenergy:endenergy) )
    822         ALLOCATE ( ddz_wall_stag(nzb_wall:nzt_wall, startenergy:endenergy) )
    823         ALLOCATE ( zw(nzb_wall:nzt_wall, startenergy:endenergy) )
    824 
    825 !--     allocate wall and roof temperature arrays
     936        ALLOCATE ( surf_usm_h%dz_wall(nzb_wall:nzt_wall+1,1:surf_usm_h%ns)     )
     937        ALLOCATE ( surf_usm_h%ddz_wall(nzb_wall:nzt_wall+1,1:surf_usm_h%ns)    )
     938        ALLOCATE ( surf_usm_h%dz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns)  )
     939        ALLOCATE ( surf_usm_h%ddz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
     940        ALLOCATE ( surf_usm_h%zw(nzb_wall:nzt_wall,1:surf_usm_h%ns)            )
     941!
     942!--     For vertical surfaces.
     943        DO  l = 0, 3
     944           ALLOCATE ( surf_usm_v(l)%dz_wall(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns)     )
     945           ALLOCATE ( surf_usm_v(l)%ddz_wall(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns)    )
     946           ALLOCATE ( surf_usm_v(l)%dz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns)  )
     947           ALLOCATE ( surf_usm_v(l)%ddz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
     948           ALLOCATE ( surf_usm_v(l)%zw(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns)            )
     949        ENDDO
     950
     951!--     allocate wall and roof temperature arrays, for horizontal walls
    826952#if defined( __nopointer )
    827         ALLOCATE ( t_surf(startenergy:endenergy) )
    828         ALLOCATE ( t_surf_p(startenergy:endenergy) )
    829         ALLOCATE ( t_wall(nzb_wall:nzt_wall+1,startenergy:endenergy) )
    830         ALLOCATE ( t_wall_p(nzb_wall:nzt_wall+1,startenergy:endenergy) )
     953        ALLOCATE ( t_surf_h(1:surf_usm_h%ns) )
     954        ALLOCATE ( t_surf_h_p(1:surf_usm_h%ns) )
     955        ALLOCATE ( t_wall_h(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
     956        ALLOCATE ( t_wall_h_p(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
     957
     958        ALLOCATE ( t_surf_h(1:surf_usm_h%ns) )
     959        ALLOCATE ( t_surf_h_p(1:surf_usm_h%ns) )
     960        ALLOCATE ( t_wall_h(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
     961        ALLOCATE ( t_wall_h_p(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
    831962#else
    832         ALLOCATE ( t_surf_1(startenergy:endenergy) )
    833         ALLOCATE ( t_surf_2(startenergy:endenergy) )
    834         ALLOCATE ( t_wall_1(nzb_wall:nzt_wall+1,startenergy:endenergy) )
    835         ALLOCATE ( t_wall_2(nzb_wall:nzt_wall+1,startenergy:endenergy) )
     963        ALLOCATE ( t_surf_h_1(1:surf_usm_h%ns) )
     964        ALLOCATE ( t_surf_h_2(1:surf_usm_h%ns) )
     965        ALLOCATE ( t_wall_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
     966        ALLOCATE ( t_wall_h_2(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
    836967
    837968!--     initial assignment of the pointers
    838         t_wall    => t_wall_1;    t_wall_p    => t_wall_2
    839         t_surf => t_surf_1; t_surf_p => t_surf_2
     969        t_wall_h    => t_wall_h_1;    t_wall_h_p    => t_wall_h_2
     970        t_surf_h => t_surf_h_1; t_surf_h_p => t_surf_h_2
    840971#endif
    841972
    842 !--     allocate intermediate timestep arrays
    843         ALLOCATE ( tt_surface_m(startenergy:endenergy) )
    844         ALLOCATE ( tt_wall_m(nzb_wall:nzt_wall+1,startenergy:endenergy) )
    845 
    846 !--     allocate wall heat flux output array
    847         ALLOCATE ( wshf(startwall:endwall) )
    848         ALLOCATE ( wshf_eb(startenergy:endenergy) )
    849         ALLOCATE ( wghf_eb(startenergy:endenergy) )
    850 
    851 !--     set inital values for prognostic quantities
    852         tt_surface_m = 0.0_wp
    853         tt_wall_m    = 0.0_wp
    854 
    855         wshf = 0.0_wp
    856         wshf_eb = 0.0_wp
    857         wghf_eb = 0.0_wp
     973!--     allocate wall and roof temperature arrays, for vertical walls
     974#if defined( __nopointer )
     975        DO  l = 0, 3
     976           ALLOCATE ( t_surf_v(l)%t(1:surf_usm_v(l)%ns) )
     977           ALLOCATE ( t_surf_v(l)%t_p(1:surf_usm_v(l)%ns) )
     978           ALLOCATE ( t_wall_v(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
     979           ALLOCATE ( t_wall_v(l)%t_p(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
     980        ENDDO
     981#else
     982        DO  l = 0, 3
     983           ALLOCATE ( t_surf_v_1(l)%t(1:surf_usm_v(l)%ns) )
     984           ALLOCATE ( t_surf_v_2(l)%t(1:surf_usm_v(l)%ns) )
     985           ALLOCATE ( t_wall_v_1(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
     986           ALLOCATE ( t_wall_v_2(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
     987        ENDDO
     988!
     989!--     initial assignment of the pointers
     990        t_wall_v    => t_wall_v_1;    t_wall_v_p    => t_wall_v_2
     991        t_surf_v => t_surf_v_1; t_surf_v_p => t_surf_v_2
     992#endif
     993!
     994!--     Allocate intermediate timestep arrays. For horizontal surfaces.
     995        ALLOCATE ( surf_usm_h%tt_surface_m(1:surf_usm_h%ns)                  )
     996        ALLOCATE ( surf_usm_h%tt_wall_m(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
     997!
     998!--     Set inital values for prognostic quantities
     999        IF ( ALLOCATED( surf_usm_h%tt_surface_m ) )  surf_usm_h%tt_surface_m = 0.0_wp
     1000        IF ( ALLOCATED( surf_usm_h%tt_wall_m    ) )  surf_usm_h%tt_wall_m    = 0.0_wp
     1001!
     1002!--     Now, for vertical surfaces
     1003        DO  l = 0, 3
     1004           ALLOCATE ( surf_usm_v(l)%tt_surface_m(1:surf_usm_v(l)%ns)                  )
     1005           ALLOCATE ( surf_usm_v(l)%tt_wall_m(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
     1006           IF ( ALLOCATED( surf_usm_v(l)%tt_surface_m ) )  surf_usm_v(l)%tt_surface_m = 0.0_wp
     1007           IF ( ALLOCATED( surf_usm_v(l)%tt_wall_m    ) )  surf_usm_v(l)%tt_wall_m    = 0.0_wp
     1008        ENDDO
     1009
     1010!--     allocate wall heat flux output array and set initial values. For horizontal surfaces
     1011!         ALLOCATE ( surf_usm_h%wshf(1:surf_usm_h%ns)    )  !can be removed
     1012        ALLOCATE ( surf_usm_h%wshf_eb(1:surf_usm_h%ns) )
     1013        ALLOCATE ( surf_usm_h%wghf_eb(1:surf_usm_h%ns) )
     1014        IF ( ALLOCATED( surf_usm_h%wshf    ) )  surf_usm_h%wshf    = 0.0_wp
     1015        IF ( ALLOCATED( surf_usm_h%wshf_eb ) )  surf_usm_h%wshf_eb = 0.0_wp
     1016        IF ( ALLOCATED( surf_usm_h%wghf_eb ) )  surf_usm_h%wghf_eb = 0.0_wp
     1017!
     1018!--     Now, for vertical surfaces
     1019        DO  l = 0, 3
     1020!            ALLOCATE ( surf_usm_v(l)%wshf(1:surf_usm_v(l)%ns)    )    ! can be removed
     1021           ALLOCATE ( surf_usm_v(l)%wshf_eb(1:surf_usm_v(l)%ns) )
     1022           ALLOCATE ( surf_usm_v(l)%wghf_eb(1:surf_usm_v(l)%ns) )
     1023           IF ( ALLOCATED( surf_usm_v(l)%wshf    ) )  surf_usm_v(l)%wshf    = 0.0_wp
     1024           IF ( ALLOCATED( surf_usm_v(l)%wshf_eb ) )  surf_usm_v(l)%wshf_eb = 0.0_wp
     1025           IF ( ALLOCATED( surf_usm_v(l)%wghf_eb ) )  surf_usm_v(l)%wghf_eb = 0.0_wp
     1026        ENDDO
    8581027       
    8591028    END SUBROUTINE usm_allocate_urban_surface
     
    8741043        CHARACTER (len=*), INTENT(IN) :: variable
    8751044 
    876         INTEGER(iwp)                                       :: i, j, k, l, ids, iwl,istat
     1045        INTEGER(iwp)                                       :: i, j, k, l, m, ids, iwl,istat
    8771046        CHARACTER (len=varnamelength)                      :: var, surfid
    8781047        INTEGER(iwp), PARAMETER                            :: nd = 5
     
    9101079                CASE ( 'usm_rad_net' )
    9111080!--                 array of complete radiation balance
    912                     IF ( .NOT.  ALLOCATED(rad_net_av) )  THEN
    913                         ALLOCATE( rad_net_av(startenergy:endenergy) )
    914                         rad_net_av = 0.0_wp
     1081                    IF ( .NOT.  ALLOCATED(surf_usm_h%rad_net_av) )  THEN
     1082                        ALLOCATE( surf_usm_h%rad_net_av(1:surf_usm_h%ns) )
     1083                        surf_usm_h%rad_net_av = 0.0_wp
    9151084                    ENDIF
     1085                    DO  l = 0, 3
     1086                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%rad_net_av) )  THEN
     1087                           ALLOCATE( surf_usm_v(l)%rad_net_av(1:surf_usm_v(l)%ns) )
     1088                           surf_usm_v(l)%rad_net_av = 0.0_wp
     1089                       ENDIF
     1090                    ENDDO
    9161091                   
    9171092                CASE ( 'usm_rad_insw' )
    9181093!--                 array of sw radiation falling to surface after i-th reflection
    919                     IF ( .NOT.  ALLOCATED(surfinsw_av) )  THEN
    920                         ALLOCATE( surfinsw_av(startenergy:endenergy) )
    921                         surfinsw_av = 0.0_wp
     1094                    IF ( .NOT.  ALLOCATED(surf_usm_h%surfinsw_av) )  THEN
     1095                        ALLOCATE( surf_usm_h%surfinsw_av(1:surf_usm_h%ns) )
     1096                        surf_usm_h%surfinsw_av = 0.0_wp
    9221097                    ENDIF
     1098                    DO  l = 0, 3
     1099                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%surfinsw_av) )  THEN
     1100                           ALLOCATE( surf_usm_v(l)%surfinsw_av(1:surf_usm_v(l)%ns) )
     1101                           surf_usm_v(l)%surfinsw_av = 0.0_wp
     1102                       ENDIF
     1103                    ENDDO
    9231104                                   
    9241105                CASE ( 'usm_rad_inlw' )
    9251106!--                 array of lw radiation falling to surface after i-th reflection
    926                     IF ( .NOT.  ALLOCATED(surfinlw_av) )  THEN
    927                         ALLOCATE( surfinlw_av(startenergy:endenergy) )
    928                         surfinlw_av = 0.0_wp
     1107                    IF ( .NOT.  ALLOCATED(surf_usm_h%surfinlw_av) )  THEN
     1108                        ALLOCATE( surf_usm_h%surfinlw_av(1:surf_usm_h%ns) )
     1109                        surf_usm_h%surfinlw_av = 0.0_wp
    9291110                    ENDIF
     1111                    DO  l = 0, 3
     1112                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%surfinlw_av) )  THEN
     1113                           ALLOCATE( surf_usm_v(l)%surfinlw_av(1:surf_usm_v(l)%ns) )
     1114                           surf_usm_v(l)%surfinlw_av = 0.0_wp
     1115                       ENDIF
     1116                    ENDDO
    9301117
    9311118                CASE ( 'usm_rad_inswdir' )
     
    9521139                CASE ( 'usm_rad_inlwdif' )
    9531140!--                 array of sw radiation falling to surface after i-th reflection
    954                     IF ( .NOT.  ALLOCATED(surfinlwdif_av) )  THEN
     1141                   IF ( .NOT.  ALLOCATED(surfinlwdif_av) )  THEN
    9551142                        ALLOCATE( surfinlwdif_av(startenergy:endenergy) )
    9561143                        surfinlwdif_av = 0.0_wp
     
    9771164                        surfoutlw_av = 0.0_wp
    9781165                    ENDIF
    979 
    9801166                CASE ( 'usm_rad_ressw' )
    9811167!--                 array of residua of sw radiation absorbed in surface after last reflection
     
    9841170                        surfins_av = 0.0_wp
    9851171                    ENDIF
    986                                     
     1172                                   
    9871173                CASE ( 'usm_rad_reslw' )
    9881174!--                 array of residua of lw radiation absorbed in surface after last reflection
     
    9941180                CASE ( 'usm_rad_hf' )
    9951181!--                 array of heat flux from radiation for surfaces after i-th reflection
    996                     IF ( .NOT.  ALLOCATED(surfhf_av) )  THEN
    997                         ALLOCATE( surfhf_av(startenergy:endenergy) )
    998                         surfhf_av = 0.0_wp
     1182                    IF ( .NOT.  ALLOCATED(surf_usm_h%surfhf_av) )  THEN
     1183                        ALLOCATE( surf_usm_h%surfhf_av(1:surf_usm_h%ns) )
     1184                        surf_usm_h%surfhf_av = 0.0_wp
    9991185                    ENDIF
     1186                    DO  l = 0, 3
     1187                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%surfhf_av) )  THEN
     1188                           ALLOCATE( surf_usm_v(l)%surfhf_av(1:surf_usm_v(l)%ns) )
     1189                           surf_usm_v(l)%surfhf_av = 0.0_wp
     1190                       ENDIF
     1191                    ENDDO
    10001192
    10011193                CASE ( 'usm_wshf' )
    10021194!--                 array of sensible heat flux from surfaces
    10031195!--                 land surfaces
    1004                     IF ( .NOT.  ALLOCATED(wshf_eb_av) )  THEN
    1005                         ALLOCATE( wshf_eb_av(startenergy:endenergy) )
    1006                         wshf_eb_av = 0.0_wp
     1196                    IF ( .NOT.  ALLOCATED(surf_usm_h%wshf_eb_av) )  THEN
     1197                        ALLOCATE( surf_usm_h%wshf_eb_av(1:surf_usm_h%ns) )
     1198                        surf_usm_h%wshf_eb_av = 0.0_wp
    10071199                    ENDIF
     1200                    DO  l = 0, 3
     1201                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%wshf_eb_av) )  THEN
     1202                           ALLOCATE( surf_usm_v(l)%wshf_eb_av(1:surf_usm_v(l)%ns) )
     1203                           surf_usm_v(l)%wshf_eb_av = 0.0_wp
     1204                       ENDIF
     1205                    ENDDO
    10081206
    10091207                CASE ( 'usm_wghf' )
    10101208!--                 array of heat flux from ground (wall, roof, land)
    1011                     IF ( .NOT.  ALLOCATED(wghf_eb_av) )  THEN
    1012                         ALLOCATE( wghf_eb_av(startenergy:endenergy) )
    1013                         wghf_eb_av = 0.0_wp
     1209                    IF ( .NOT.  ALLOCATED(surf_usm_h%wghf_eb_av) )  THEN
     1210                        ALLOCATE( surf_usm_h%wghf_eb_av(1:surf_usm_h%ns) )
     1211                        surf_usm_h%wghf_eb_av = 0.0_wp
    10141212                    ENDIF
     1213                    DO  l = 0, 3
     1214                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%wghf_eb_av) )  THEN
     1215                           ALLOCATE( surf_usm_v(l)%wghf_eb_av(1:surf_usm_v(l)%ns) )
     1216                           surf_usm_v(l)%wghf_eb_av = 0.0_wp
     1217                       ENDIF
     1218                    ENDDO
    10151219
    10161220                CASE ( 'usm_t_surf' )
    10171221!--                 surface temperature for surfaces
    1018                     IF ( .NOT.  ALLOCATED(t_surf_av) )  THEN
    1019                         ALLOCATE( t_surf_av(startenergy:endenergy) )
    1020                         t_surf_av = 0.0_wp
     1222                    IF ( .NOT.  ALLOCATED(surf_usm_h%t_surf_av) )  THEN
     1223                        ALLOCATE( surf_usm_h%t_surf_av(1:surf_usm_h%ns) )
     1224                        surf_usm_h%t_surf_av = 0.0_wp
    10211225                    ENDIF
     1226                    DO  l = 0, 3
     1227                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_surf_av) )  THEN
     1228                           ALLOCATE( surf_usm_v(l)%t_surf_av(1:surf_usm_v(l)%ns) )
     1229                           surf_usm_v(l)%t_surf_av = 0.0_wp
     1230                       ENDIF
     1231                    ENDDO
    10221232
    10231233                CASE ( 'usm_t_wall' )
    10241234!--                 wall temperature for iwl layer of walls and land
    1025                     IF ( .NOT.  ALLOCATED(t_wall_av) )  THEN
    1026                         ALLOCATE( t_wall_av(nzb_wall:nzt_wall,startenergy:endenergy) )
    1027                         t_wall_av = 0.0_wp
     1235                    IF ( .NOT.  ALLOCATED(surf_usm_h%t_wall_av) )  THEN
     1236                        ALLOCATE( surf_usm_h%t_wall_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
     1237                        surf_usm_h%t_wall_av = 0.0_wp
    10281238                    ENDIF
     1239                    DO  l = 0, 3
     1240                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_wall_av) )  THEN
     1241                           ALLOCATE( surf_usm_v(l)%t_wall_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
     1242                           surf_usm_v(l)%t_wall_av = 0.0_wp
     1243                       ENDIF
     1244                    ENDDO
    10291245
    10301246               CASE DEFAULT
     
    10391255                CASE ( 'usm_rad_net' )
    10401256!--                 array of complete radiation balance
    1041                     DO l = startenergy, endenergy
    1042                         IF ( surfl(id,l) == ids )  THEN
    1043                             rad_net_av(l) = rad_net_av(l) + rad_net_l(l)
    1044                         ENDIF
     1257                    DO  m = 1, surf_usm_h%ns
     1258                       surf_usm_h%rad_net_av(m) =                              &
     1259                                          surf_usm_h%rad_net_av(m) +           &
     1260                                          surf_usm_h%rad_net_l(m)
     1261                    ENDDO
     1262                    DO  l = 0, 3
     1263                       DO  m = 1, surf_usm_v(l)%ns
     1264                          surf_usm_v(l)%rad_net_av(m) =                        &
     1265                                          surf_usm_v(l)%rad_net_av(m) +        &
     1266                                          surf_usm_v(l)%rad_net_l(m)
     1267                       ENDDO
    10451268                    ENDDO
    10461269                   
     
    10851308                        ENDIF
    10861309                    ENDDO
     1310
    10871311                   
    10881312                CASE ( 'usm_rad_inlwdif' )
     
    10901314                    DO l = startenergy, endenergy
    10911315                        IF ( surfl(id,l) == ids )  THEN
     1316                            surfinswref_av(l) = surfinswref_av(l) + surfinsw(l) - &
     1317                                                surfinswdir(l) - surfinswdif(l)
     1318                        ENDIF
     1319                    ENDDO
     1320!                     
     1321                CASE ( 'usm_rad_inlwref' )
     1322!--                 array of lw radiation falling to surface from reflections
     1323                    DO l = startenergy, endenergy
     1324                        IF ( surfl(id,l) == ids )  THEN
    10921325                            surfinlwdif_av(l) = surfinlwdif_av(l) + surfinlwdif(l)
    10931326                        ENDIF
    10941327                    ENDDO
    10951328                   
    1096                 CASE ( 'usm_rad_inlwref' )
    1097 !--                 array of lw radiation falling to surface from reflections
     1329                CASE ( 'usm_rad_outsw' )
     1330!--                 array of sw radiation emitted from surface after i-th reflection
    10981331                    DO l = startenergy, endenergy
    10991332                        IF ( surfl(id,l) == ids )  THEN
     
    11031336                    ENDDO
    11041337                   
    1105                 CASE ( 'usm_rad_outsw' )
    1106 !--                 array of sw radiation emitted from surface after i-th reflection
    1107                     DO l = startenergy, endenergy
    1108                         IF ( surfl(id,l) == ids )  THEN
    1109                             surfoutsw_av(l) = surfoutsw_av(l) + surfoutsw(l)
    1110                         ENDIF
    1111                     ENDDO
    1112                    
    11131338                CASE ( 'usm_rad_outlw' )
    11141339!--                 array of lw radiation emitted from surface after i-th reflection
    11151340                    DO l = startenergy, endenergy
    11161341                        IF ( surfl(id,l) == ids )  THEN
    1117                             surfoutlw_av(l) = surfoutlw_av(l) + surfoutlw(l)
     1342                            surfoutsw_av(l) = surfoutsw_av(l) + surfoutsw(l)
    11181343                        ENDIF
    11191344                    ENDDO
     
    11231348                    DO l = startenergy, endenergy
    11241349                        IF ( surfl(id,l) == ids )  THEN
    1125                             surfins_av(l) = surfins_av(l) + surfins(l)
     1350                            surfoutlw_av(l) = surfoutlw_av(l) + surfoutlw(l)
    11261351                        ENDIF
    11271352                    ENDDO
     
    11311356                    DO l = startenergy, endenergy
    11321357                        IF ( surfl(id,l) == ids )  THEN
    1133                             surfinl_av(l) = surfinl_av(l) + surfinl(l)
     1358                            surfins_av(l) = surfins_av(l) + surfins(l)
    11341359                        ENDIF
    11351360                    ENDDO
     
    11371362                CASE ( 'usm_rad_hf' )
    11381363!--                 array of heat flux from radiation for surfaces after i-th reflection
    1139                     DO l = startenergy, endenergy
    1140                         IF ( surfl(id,l) == ids )  THEN
    1141                             surfhf_av(l) = surfhf_av(l) + surfhf(l)
    1142                         ENDIF
     1364                    DO  m = 1, surf_usm_h%ns
     1365                       surf_usm_h%surfhf_av(m) =                               &
     1366                                          surf_usm_h%surfhf_av(m) +            &
     1367                                          surf_usm_h%surfhf(m)
     1368                    ENDDO
     1369                    DO  l = 0, 3
     1370                       DO  m = 1, surf_usm_v(l)%ns
     1371                          surf_usm_v(l)%surfhf_av(m) =                         &
     1372                                          surf_usm_v(l)%surfhf_av(m) +         &
     1373                                          surf_usm_v(l)%surfhf(m)
     1374                       ENDDO
    11431375                    ENDDO
    11441376                   
    11451377                CASE ( 'usm_wshf' )
    11461378!--                 array of sensible heat flux from surfaces (land, roof, wall)
    1147                     DO l = startenergy, endenergy
    1148                         IF ( surfl(id,l) == ids )  THEN
    1149                             wshf_eb_av(l) = wshf_eb_av(l) + wshf_eb(l)
    1150                         ENDIF
     1379                    DO  m = 1, surf_usm_h%ns
     1380                       surf_usm_h%wshf_eb_av(m) =                              &
     1381                                          surf_usm_h%wshf_eb_av(m) +           &
     1382                                          surf_usm_h%wshf_eb(m)
     1383                    ENDDO
     1384                    DO  l = 0, 3
     1385                       DO  m = 1, surf_usm_v(l)%ns
     1386                          surf_usm_v(l)%wshf_eb_av(m) =                        &
     1387                                          surf_usm_v(l)%wshf_eb_av(m) +        &
     1388                                          surf_usm_v(l)%wshf_eb(m)
     1389                       ENDDO
    11511390                    ENDDO
    11521391                   
    11531392                CASE ( 'usm_wghf' )
    11541393!--                 array of heat flux from ground (wall, roof, land)
    1155                     DO l = startenergy, endenergy
    1156                         IF ( surfl(id,l) == ids )  THEN
    1157                             wghf_eb_av(l) = wghf_eb_av(l) + wghf_eb(l)
    1158                         ENDIF
     1394                    DO  m = 1, surf_usm_h%ns
     1395                       surf_usm_h%wghf_eb_av(m) =                              &
     1396                                          surf_usm_h%wghf_eb_av(m) +           &
     1397                                          surf_usm_h%wghf_eb(m)
     1398                    ENDDO
     1399                    DO  l = 0, 3
     1400                       DO  m = 1, surf_usm_v(l)%ns
     1401                          surf_usm_v(l)%wghf_eb_av(m) =                        &
     1402                                          surf_usm_v(l)%wghf_eb_av(m) +        &
     1403                                          surf_usm_v(l)%wghf_eb(m)
     1404                       ENDDO
    11591405                    ENDDO
    11601406                   
    11611407                CASE ( 'usm_t_surf' )
    11621408!--                 surface temperature for surfaces
    1163                     DO l = startenergy, endenergy
    1164                         IF ( surfl(id,l) == ids )  THEN
    1165                             t_surf_av(l) = t_surf_av(l) + t_surf(l)
    1166                         ENDIF
     1409                    DO  m = 1, surf_usm_h%ns
     1410                       surf_usm_h%t_surf_av(m) =                               &
     1411                                          surf_usm_h%t_surf_av(m) +            &
     1412                                          t_surf_h(m)
     1413                    ENDDO
     1414                    DO  l = 0, 3
     1415                       DO  m = 1, surf_usm_v(l)%ns
     1416                          surf_usm_v(l)%t_surf_av(m) =                         &
     1417                                          surf_usm_v(l)%t_surf_av(m) +         &
     1418                                          t_surf_v(l)%t(m)
     1419                       ENDDO
    11671420                    ENDDO
    11681421                   
    11691422                CASE ( 'usm_t_wall' )
    11701423!--                 wall temperature for  iwl layer of walls and land
    1171                     DO l = startenergy, endenergy
    1172                         IF ( surfl(id,l) == ids )  THEN
    1173                             t_wall_av(iwl, l) = t_wall_av(iwl,l) + t_wall(iwl, l)
    1174                         ENDIF
     1424                    DO  m = 1, surf_usm_h%ns
     1425                       surf_usm_h%t_wall_av(iwl,m) =                           &
     1426                                          surf_usm_h%t_wall_av(iwl,m) +        &
     1427                                          t_wall_h(iwl,m)
     1428                    ENDDO
     1429                    DO  l = 0, 3
     1430                       DO  m = 1, surf_usm_v(l)%ns
     1431                          surf_usm_v(l)%t_wall_av(iwl,m) =                     &
     1432                                          surf_usm_v(l)%t_wall_av(iwl,m) +     &
     1433                                          t_wall_v(l)%t(iwl,m)
     1434                       ENDDO
    11751435                    ENDDO
    11761436                   
     
    11861446                CASE ( 'usm_rad_net' )
    11871447!--                 array of complete radiation balance
    1188                     DO l = startenergy, endenergy
    1189                         IF ( surfl(id,l) == ids )  THEN
    1190                             rad_net_av(l) = rad_net_av(l) / REAL( average_count_3d, kind=wp )
    1191                         ENDIF
     1448                    DO  m = 1, surf_usm_h%ns
     1449                       surf_usm_h%rad_net_av(m) =                              &
     1450                                          surf_usm_h%rad_net_av(m) /           &
     1451                                          REAL( average_count_3d, kind=wp )
     1452                    ENDDO
     1453                    DO  l = 0, 3
     1454                       DO  m = 1, surf_usm_v(l)%ns
     1455                          surf_usm_v(l)%rad_net_av(m) =                        &
     1456                                          surf_usm_v(l)%rad_net_av(m) /        &
     1457                                          REAL( average_count_3d, kind=wp )
     1458                       ENDDO
    11921459                    ENDDO
    11931460                   
     
    11991466                        ENDIF
    12001467                    ENDDO
    1201                                     
     1468                             
    12021469                CASE ( 'usm_rad_inlw' )
    12031470!--                 array of lw radiation falling to surface after i-th reflection
     
    12071474                        ENDIF
    12081475                    ENDDO
    1209 
     1476                   
    12101477                CASE ( 'usm_rad_inswdir' )
    12111478!--                 array of direct sw radiation falling to surface from sun
     
    12151482                        ENDIF
    12161483                    ENDDO
    1217 
     1484                   
    12181485                CASE ( 'usm_rad_inswdif' )
    12191486!--                 array of difusion sw radiation falling to surface from sky and borders of the domain
     
    12231490                        ENDIF
    12241491                    ENDDO
    1225 
     1492                   
    12261493                CASE ( 'usm_rad_inswref' )
    12271494!--                 array of sw radiation falling to surface from reflections
     
    12311498                        ENDIF
    12321499                    ENDDO
    1233 
     1500                   
    12341501                CASE ( 'usm_rad_inlwdif' )
    12351502!--                 array of sw radiation falling to surface after i-th reflection
     
    12391506                        ENDIF
    12401507                    ENDDO
    1241 
     1508                   
    12421509                CASE ( 'usm_rad_inlwref' )
    12431510!--                 array of lw radiation falling to surface from reflections
     
    12471514                        ENDIF
    12481515                    ENDDO
    1249 
     1516                   
    12501517                CASE ( 'usm_rad_outsw' )
    12511518!--                 array of sw radiation emitted from surface after i-th reflection
     
    12551522                        ENDIF
    12561523                    ENDDO
    1257 
     1524                   
    12581525                CASE ( 'usm_rad_outlw' )
    12591526!--                 array of lw radiation emitted from surface after i-th reflection
     
    12631530                        ENDIF
    12641531                    ENDDO
    1265 
     1532                   
    12661533                CASE ( 'usm_rad_ressw' )
    12671534!--                 array of residua of sw radiation absorbed in surface after last reflection
     
    12821549                CASE ( 'usm_rad_hf' )
    12831550!--                 array of heat flux from radiation for surfaces after i-th reflection
    1284                     DO l = startenergy, endenergy
    1285                         IF ( surfl(id,l) == ids )  THEN
    1286                             surfhf_av(l) = surfhf_av(l) / REAL( average_count_3d, kind=wp )
    1287                         ENDIF
    1288                     ENDDO
    1289 
     1551                    DO  m = 1, surf_usm_h%ns
     1552                       surf_usm_h%surfhf_av(m) =                               &
     1553                                          surf_usm_h%surfhf_av(m) /            &
     1554                                          REAL( average_count_3d, kind=wp )
     1555                    ENDDO
     1556                    DO  l = 0, 3
     1557                       DO  m = 1, surf_usm_v(l)%ns
     1558                          surf_usm_v(l)%surfhf_av(m) =                         &
     1559                                          surf_usm_v(l)%surfhf_av(m) /         &
     1560                                          REAL( average_count_3d, kind=wp )
     1561                       ENDDO
     1562                    ENDDO
     1563                   
    12901564                CASE ( 'usm_wshf' )
    12911565!--                 array of sensible heat flux from surfaces (land, roof, wall)
    1292                     DO l = startenergy, endenergy
    1293                         IF ( surfl(id,l) == ids )  THEN
    1294                             wshf_eb_av(l) = wshf_eb_av(l) / REAL( average_count_3d, kind=wp )
    1295                         ENDIF
    1296                     ENDDO
    1297 
     1566                    DO  m = 1, surf_usm_h%ns
     1567                       surf_usm_h%wshf_eb_av(m) =                              &
     1568                                          surf_usm_h%wshf_eb_av(m) /           &
     1569                                          REAL( average_count_3d, kind=wp )
     1570                    ENDDO
     1571                    DO  l = 0, 3
     1572                       DO  m = 1, surf_usm_v(l)%ns
     1573                          surf_usm_v(l)%wshf_eb_av(m) =                        &
     1574                                          surf_usm_v(l)%wshf_eb_av(m) /        &
     1575                                          REAL( average_count_3d, kind=wp )
     1576                       ENDDO
     1577                    ENDDO
     1578                   
    12981579                CASE ( 'usm_wghf' )
    12991580!--                 array of heat flux from ground (wall, roof, land)
    1300                     DO l = startenergy, endenergy
    1301                         IF ( surfl(id,l) == ids )  THEN
    1302                             wghf_eb_av(l) = wghf_eb_av(l) / REAL( average_count_3d, kind=wp )
    1303                         ENDIF
    1304                     ENDDO
    1305 
     1581                    DO  m = 1, surf_usm_h%ns
     1582                       surf_usm_h%wghf_eb_av(m) =                              &
     1583                                          surf_usm_h%wghf_eb_av(m) /           &
     1584                                          REAL( average_count_3d, kind=wp )
     1585                    ENDDO
     1586                    DO  l = 0, 3
     1587                       DO  m = 1, surf_usm_v(l)%ns
     1588                          surf_usm_v(l)%wghf_eb_av(m) =                        &
     1589                                          surf_usm_v(l)%wghf_eb_av(m) /        &
     1590                                          REAL( average_count_3d, kind=wp )
     1591                       ENDDO
     1592                    ENDDO
     1593                   
    13061594                CASE ( 'usm_t_surf' )
    13071595!--                 surface temperature for surfaces
    1308                     DO l = startenergy, endenergy
    1309                         IF ( surfl(id,l) == ids )  THEN
    1310                             t_surf_av(l) = t_surf_av(l) / REAL( average_count_3d, kind=wp )
    1311                         ENDIF
    1312                     ENDDO
    1313 
     1596                    DO  m = 1, surf_usm_h%ns
     1597                       surf_usm_h%t_surf_av(m) =                               &
     1598                                          surf_usm_h%t_surf_av(m) /            &
     1599                                          REAL( average_count_3d, kind=wp )
     1600                    ENDDO
     1601                    DO  l = 0, 3
     1602                       DO  m = 1, surf_usm_v(l)%ns
     1603                          surf_usm_v(l)%t_surf_av(m) =                         &
     1604                                          surf_usm_v(l)%t_surf_av(m) /         &
     1605                                          REAL( average_count_3d, kind=wp )
     1606                       ENDDO
     1607                    ENDDO
     1608                   
    13141609                CASE ( 'usm_t_wall' )
    13151610!--                 wall temperature for  iwl layer of walls and land
    1316                     DO l = startenergy, endenergy
    1317                         IF ( surfl(id,l) == ids )  THEN
    1318                             t_wall_av(iwl, l) = t_wall_av(iwl,l) / REAL( average_count_3d, kind=wp )
    1319                         ENDIF
     1611                    DO  m = 1, surf_usm_h%ns
     1612                       surf_usm_h%t_wall_av(iwl,m) =                           &
     1613                                          surf_usm_h%t_wall_av(iwl,m) /        &
     1614                                          REAL( average_count_3d, kind=wp )
     1615                    ENDDO
     1616                    DO  l = 0, 3
     1617                       DO  m = 1, surf_usm_v(l)%ns
     1618                          surf_usm_v(l)%t_wall_av(iwl,m) =                     &
     1619                                          surf_usm_v(l)%t_wall_av(iwl,m) /     &
     1620                                          REAL( average_count_3d, kind=wp )
     1621                       ENDDO
    13201622                    ENDDO
    13211623
     
    14931795        TYPE(c_ptr)                                 :: lad_s_rma_p     !< allocated c pointer
    14941796        INTEGER(kind=MPI_ADDRESS_KIND)              :: size_lad_rma
    1495    
     1797!   
    14961798!--     calculation of the SVF
    14971799        CALL location_message( '    calculation of SVF and CSF', .TRUE. )
    1498 
    1499 
     1800!
    15001801!--     precalculate face areas for different face directions using normal vector
    15011802        DO d = 0, 9
     
    15241825#if defined( __parallel )
    15251826        ALLOCATE( nzterrl(nys:nyn,nxl:nxr) )
    1526         nzterrl = nzb_s_inner(nys:nyn,nxl:nxr)
     1827        nzterrl = MAXLOC(                                                      &
     1828                          MERGE( 1, 0,                                         &
     1829                                 BTEST( wall_flags_0(:,nys:nyn,nxl:nxr), 12 )  &
     1830                               ), DIM = 1                                      &
     1831                        ) - 1  ! = nzb_s_inner(nys:nyn,nxl:nxr)
    15271832        CALL MPI_AllGather( nzterrl, nnx*nny, MPI_INTEGER, &
    15281833                            nzterr, nnx*nny, MPI_INTEGER, comm2d, ierr )
    15291834        DEALLOCATE(nzterrl)
    15301835#else
    1531         nzterr = RESHAPE( nzb_s_inner(nys:nyn,nxl:nxr), (/(nx+1)*(ny+1)/) )
     1836        nzterr = RESHAPE( MAXLOC(                                              &
     1837                          MERGE( 1, 0,                                         &
     1838                                 BTEST( wall_flags_0(:,nys:nyn,nxl:nxr), 12 )  &
     1839                               ), DIM = 1                                      &
     1840                                ) - 1,                                         &
     1841                          (/(nx+1)*(ny+1)/)                                    &
     1842                        )
    15321843#endif
    15331844        IF ( plant_canopy )  THEN
     
    15791890            DO i = nxl, nxr
    15801891                DO j = nys, nyn
    1581                     k = nzb_s_inner(j, i)
     1892                    k = MAXLOC(                                                &
     1893                                MERGE( 1, 0,                                   &
     1894                                       BTEST( wall_flags_0(:,j,i), 12 )        &
     1895                                     ), DIM = 1                                &
     1896                              ) - 1
     1897
    15821898                    usm_lad(k:nzut, j, i) = lad_s(0:nzut-k, j, i)
    15831899                ENDDO
     
    19722288        CALL message( 'init_urban_surface', 'PA0502', 2, 2, 0, 6, 0 )
    19732289       
    1974 
    19752290    END SUBROUTINE usm_calc_svf
    19762291
     
    20892404        INTEGER(iwp), DIMENSION(0:nd-1)                        :: dirend
    20902405        INTEGER(iwp)                                           :: ids,isurf,isvf,isurfs,isurflt
    2091         INTEGER(iwp)                                           :: is,js,ks,i,j,k,iwl,istat
     2406        INTEGER(iwp)                                           :: is,js,ks,i,j,k,iwl,istat, l, m
     2407        INTEGER(iwp)                                           ::  k_topo    !< topography top index
    20922408
    20932409        dirstart = (/ startland, startwall, startwall, startwall, startwall /)
     
    21392455          CASE ( 'usm_surfz' )
    21402456!--           array of lw radiation falling to local surface after i-th reflection
    2141               DO isurf = dirstart(ids), dirend(ids)
    2142                   IF ( surfl(id,isurf) == ids )  THEN
    2143                       IF ( surfl(id,isurf) == iroof )  THEN
    2144                           temp_pf(0,surfl(iy,isurf),surfl(ix,isurf)) =             &
    2145                                   max(temp_pf(0,surfl(iy,isurf),surfl(ix,isurf)),  &
    2146                                       REAL(surfl(iz,isurf),wp))
    2147                       ELSE
    2148                           temp_pf(0,surfl(iy,isurf),surfl(ix,isurf)) =             &
    2149                                   max(temp_pf(0,surfl(iy,isurf),surfl(ix,isurf)),  &
    2150                                       REAL(surfl(iz,isurf),wp)+1.0_wp)
    2151                       ENDIF
    2152                   ENDIF
     2457              DO  m = 1, surf_usm_h%ns
     2458                 i = surf_usm_h%i(m)
     2459                 j = surf_usm_h%j(m)
     2460                 k = surf_usm_h%k(m)
     2461                 temp_pf(0,j,i) = MAX( temp_pf(0,j,i), REAL( k, kind=wp) )
     2462              ENDDO
     2463              DO  l = 0, 3
     2464                 DO  m = 1, surf_usm_v(l)%ns
     2465                    i = surf_usm_v(l)%i(m)
     2466                    j = surf_usm_v(l)%j(m)
     2467                    k = surf_usm_v(l)%k(m)
     2468                    temp_pf(0,j,i) = MAX( temp_pf(0,j,i), REAL( k, kind=wp) + 1.0_wp )
     2469                 ENDDO
    21532470              ENDDO
    21542471
    21552472          CASE ( 'usm_surfcat' )
    21562473!--           surface category
    2157               DO isurf = dirstart(ids), dirend(ids)
    2158                  IF ( surfl(id,isurf) == ids )  THEN
    2159                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surface_types(isurf)
    2160                  ENDIF
     2474              DO  m = 1, surf_usm_h%ns
     2475                 i = surf_usm_h%i(m)
     2476                 j = surf_usm_h%j(m)
     2477                 k = surf_usm_h%k(m)
     2478                 temp_pf(k,j,i) = surf_usm_h%surface_types(m)
     2479              ENDDO
     2480              DO  l = 0, 3
     2481                 DO  m = 1, surf_usm_v(l)%ns
     2482                    i = surf_usm_v(l)%i(m)
     2483                    j = surf_usm_v(l)%j(m)
     2484                    k = surf_usm_v(l)%k(m)
     2485                    temp_pf(k,j,i) = surf_usm_v(l)%surface_types(m)
     2486                 ENDDO
    21612487              ENDDO
    21622488             
    21632489          CASE ( 'usm_surfalb' )
    21642490!--           surface albedo
    2165               DO isurf = dirstart(ids), dirend(ids)
    2166                  IF ( surfl(id,isurf) == ids )  THEN
    2167                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = albedo_surf(isurf)
    2168                  ENDIF
     2491              DO  m = 1, surf_usm_h%ns
     2492                 i = surf_usm_h%i(m)
     2493                 j = surf_usm_h%j(m)
     2494                 k = surf_usm_h%k(m)
     2495                 temp_pf(k,j,i) = surf_usm_h%albedo_surf(m)
     2496              ENDDO
     2497              DO  l = 0, 3
     2498                 DO  m = 1, surf_usm_v(l)%ns
     2499                    i = surf_usm_v(l)%i(m)
     2500                    j = surf_usm_v(l)%j(m)
     2501                    k = surf_usm_v(l)%k(m)
     2502                    temp_pf(k,j,i) = surf_usm_v(l)%albedo_surf(m)
     2503                 ENDDO
    21692504              ENDDO
    21702505             
    21712506          CASE ( 'usm_surfemis' )
    21722507!--           surface albedo
    2173               DO isurf = dirstart(ids), dirend(ids)
    2174                  IF ( surfl(id,isurf) == ids )  THEN
    2175                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = emiss_surf(isurf)
    2176                  ENDIF
     2508              DO  m = 1, surf_usm_h%ns
     2509                 i = surf_usm_h%i(m)
     2510                 j = surf_usm_h%j(m)
     2511                 k = surf_usm_h%k(m)
     2512                 temp_pf(k,j,i) = surf_usm_h%emiss_surf(m)
    21772513              ENDDO
    2178              
     2514              DO  l = 0, 3
     2515                 DO  m = 1, surf_usm_v(l)%ns
     2516                    i = surf_usm_v(l)%i(m)
     2517                    j = surf_usm_v(l)%j(m)
     2518                    k = surf_usm_v(l)%k(m)
     2519                    temp_pf(k,j,i) = surf_usm_v(l)%emiss_surf(m)
     2520                 ENDDO
     2521              ENDDO
     2522!
     2523!-- Not adjusted so far             
    21792524          CASE ( 'usm_svf', 'usm_dif' )
    21802525!--           shape view factors or iradiance factors to selected surface
     
    21972542          CASE ( 'usm_rad_net' )
    21982543!--           array of complete radiation balance
    2199               DO isurf = dirstart(ids), dirend(ids)
    2200                  IF ( surfl(id,isurf) == ids )  THEN
    2201                    IF ( av == 0 )  THEN
    2202                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = rad_net_l(isurf)
    2203                    ELSE
    2204                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = rad_net_av(isurf)
    2205                    ENDIF
    2206                  ENDIF
    2207               ENDDO
     2544              IF ( av == 0 )  THEN
     2545                 DO  m = 1, surf_usm_h%ns
     2546                    i = surf_usm_h%i(m)
     2547                    j = surf_usm_h%j(m)
     2548                    k = surf_usm_h%k(m)
     2549                    temp_pf(k,j,i) = surf_usm_h%rad_net_l(m)
     2550                 ENDDO
     2551                 DO  l = 0, 3
     2552                    DO  m = 1, surf_usm_v(l)%ns
     2553                       i = surf_usm_v(l)%i(m)
     2554                       j = surf_usm_v(l)%j(m)
     2555                       k = surf_usm_v(l)%k(m)
     2556                       temp_pf(k,j,i) = surf_usm_v(l)%rad_net_l(m)
     2557                    ENDDO
     2558                 ENDDO
     2559              ELSE
     2560                 DO  m = 1, surf_usm_h%ns
     2561                    i = surf_usm_h%i(m)
     2562                    j = surf_usm_h%j(m)
     2563                    k = surf_usm_h%k(m)
     2564                    temp_pf(k,j,i) = surf_usm_h%rad_net_av(m)
     2565                 ENDDO
     2566                 DO  l = 0, 3
     2567                    DO  m = 1, surf_usm_v(l)%ns
     2568                       i = surf_usm_v(l)%i(m)
     2569                       j = surf_usm_v(l)%j(m)
     2570                       k = surf_usm_v(l)%k(m)
     2571                       temp_pf(k,j,i) = surf_usm_v(l)%rad_net_av(m)
     2572                    ENDDO
     2573                 ENDDO
     2574              ENDIF
    22082575
    22092576          CASE ( 'usm_rad_insw' )
     
    22682635              ENDDO
    22692636
    2270           CASE ( 'usm_rad_inlwdif' )
    2271 !--           array of sw radiation falling to surface after i-th reflection
    2272               DO isurf = dirstart(ids), dirend(ids)
    2273                  IF ( surfl(id,isurf) == ids )  THEN
    2274                    IF ( av == 0 )  THEN
    2275                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlwdif(isurf)
    2276                    ELSE
    2277                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlwdif_av(isurf)
    2278                    ENDIF
    2279                  ENDIF
    2280               ENDDO
    2281 
    22822637          CASE ( 'usm_rad_inlwref' )
    22832638!--           array of lw radiation falling to surface from reflections
     
    23392694                 ENDIF
    23402695              ENDDO
    2341 
     2696 
    23422697          CASE ( 'usm_rad_hf' )
    23432698!--           array of heat flux from radiation for surfaces after all reflections
    2344               DO isurf = dirstart(ids), dirend(ids)
    2345                  IF ( surfl(id,isurf) == ids )  THEN
    2346                    IF ( av == 0 )  THEN
    2347                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfhf(isurf)
    2348                    ELSE
    2349                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfhf_av(isurf)
    2350                    ENDIF
    2351                  ENDIF
    2352               ENDDO
    2353 
     2699              IF ( av == 0 )  THEN
     2700                 DO  m = 1, surf_usm_h%ns
     2701                    i = surf_usm_h%i(m)
     2702                    j = surf_usm_h%j(m)
     2703                    k = surf_usm_h%k(m)
     2704                    temp_pf(k,j,i) = surf_usm_h%surfhf(m)
     2705                 ENDDO
     2706                 DO  l = 0, 3
     2707                    DO  m = 1, surf_usm_v(l)%ns
     2708                       i = surf_usm_v(l)%i(m)
     2709                       j = surf_usm_v(l)%j(m)
     2710                       k = surf_usm_v(l)%k(m)
     2711                       temp_pf(k,j,i) = surf_usm_v(l)%surfhf(m)
     2712                    ENDDO
     2713                 ENDDO
     2714              ELSE
     2715                 DO  m = 1, surf_usm_h%ns
     2716                    i = surf_usm_h%i(m)
     2717                    j = surf_usm_h%j(m)
     2718                    k = surf_usm_h%k(m)
     2719                    temp_pf(k,j,i) = surf_usm_h%surfhf_av(m)
     2720                 ENDDO
     2721                 DO  l = 0, 3
     2722                    DO  m = 1, surf_usm_v(l)%ns
     2723                       i = surf_usm_v(l)%i(m)
     2724                       j = surf_usm_v(l)%j(m)
     2725                       k = surf_usm_v(l)%k(m)
     2726                       temp_pf(k,j,i) = surf_usm_v(l)%surfhf_av(m)
     2727                    ENDDO
     2728                 ENDDO
     2729              ENDIF
     2730 
    23542731          CASE ( 'usm_wshf' )
    23552732!--           array of sensible heat flux from surfaces
    2356 !--           horizontal surfaces
    2357               DO isurf = dirstart(ids), dirend(ids)
    2358                  IF ( surfl(id,isurf) == ids )  THEN
    2359                    IF ( av == 0 )  THEN
    2360                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = wshf_eb(isurf)
    2361                    ELSE
    2362                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = wshf_eb_av(isurf)
    2363                    ENDIF
    2364                  ENDIF
    2365               ENDDO
     2733              IF ( av == 0 )  THEN
     2734                 DO  m = 1, surf_usm_h%ns
     2735                    i = surf_usm_h%i(m)
     2736                    j = surf_usm_h%j(m)
     2737                    k = surf_usm_h%k(m)
     2738                    temp_pf(k,j,i) = surf_usm_h%wshf_eb(m)
     2739                 ENDDO
     2740                 DO  l = 0, 3
     2741                    DO  m = 1, surf_usm_v(l)%ns
     2742                       i = surf_usm_v(l)%i(m)
     2743                       j = surf_usm_v(l)%j(m)
     2744                       k = surf_usm_v(l)%k(m)
     2745                       temp_pf(k,j,i) = surf_usm_v(l)%wshf_eb(m)
     2746                    ENDDO
     2747                 ENDDO
     2748              ELSE
     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%wshf_eb_av(m)
     2754                 ENDDO
     2755                 DO  l = 0, 3
     2756                    DO  m = 1, surf_usm_v(l)%ns
     2757                       i = surf_usm_v(l)%i(m)
     2758                       j = surf_usm_v(l)%j(m)
     2759                       k = surf_usm_v(l)%k(m)
     2760                       temp_pf(k,j,i) = surf_usm_v(l)%wshf_eb_av(m)
     2761                    ENDDO
     2762                 ENDDO
     2763              ENDIF
     2764
    23662765
    23672766          CASE ( 'usm_wghf' )
    23682767!--           array of heat flux from ground (land, wall, roof)
    2369               DO isurf = dirstart(ids), dirend(ids)
    2370                  IF ( surfl(id,isurf) == ids )  THEN
    2371                    IF ( av == 0 )  THEN
    2372                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = wghf_eb(isurf)
    2373                    ELSE
    2374                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = wghf_eb_av(isurf)
    2375                    ENDIF
    2376                  ENDIF
    2377               ENDDO
     2768              IF ( av == 0 )  THEN
     2769                 DO  m = 1, surf_usm_h%ns
     2770                    i = surf_usm_h%i(m)
     2771                    j = surf_usm_h%j(m)
     2772                    k = surf_usm_h%k(m)
     2773                    temp_pf(k,j,i) = surf_usm_h%wghf_eb(m)
     2774                 ENDDO
     2775                 DO  l = 0, 3
     2776                    DO  m = 1, surf_usm_v(l)%ns
     2777                       i = surf_usm_v(l)%i(m)
     2778                       j = surf_usm_v(l)%j(m)
     2779                       k = surf_usm_v(l)%k(m)
     2780                       temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb(m)
     2781                    ENDDO
     2782                 ENDDO
     2783              ELSE
     2784                 DO  m = 1, surf_usm_h%ns
     2785                    i = surf_usm_h%i(m)
     2786                    j = surf_usm_h%j(m)
     2787                    k = surf_usm_h%k(m)
     2788                    temp_pf(k,j,i) = surf_usm_h%wghf_eb_av(m)
     2789                 ENDDO
     2790                 DO  l = 0, 3
     2791                    DO  m = 1, surf_usm_v(l)%ns
     2792                       i = surf_usm_v(l)%i(m)
     2793                       j = surf_usm_v(l)%j(m)
     2794                       k = surf_usm_v(l)%k(m)
     2795                       temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_av(m)
     2796                    ENDDO
     2797                 ENDDO
     2798              ENDIF
    23782799
    23792800          CASE ( 'usm_t_surf' )
    23802801!--           surface temperature for surfaces
    2381               DO isurf = max(startenergy,dirstart(ids)), min(endenergy,dirend(ids))
    2382                  IF ( surfl(id,isurf) == ids )  THEN
    2383                    IF ( av == 0 )  THEN
    2384                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = t_surf(isurf)
    2385                    ELSE
    2386                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = t_surf_av(isurf)
    2387                    ENDIF
    2388                  ENDIF
    2389               ENDDO
     2802              IF ( av == 0 )  THEN
     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) = t_surf_h(m)
     2808                 ENDDO
     2809                 DO  l = 0, 3
     2810                    DO  m = 1, surf_usm_v(l)%ns
     2811                       i = surf_usm_v(l)%i(m)
     2812                       j = surf_usm_v(l)%j(m)
     2813                       k = surf_usm_v(l)%k(m)
     2814                       temp_pf(k,j,i) = t_surf_v(l)%t(m)
     2815                    ENDDO
     2816                 ENDDO
     2817              ELSE
     2818                 DO  m = 1, surf_usm_h%ns
     2819                    i = surf_usm_h%i(m)
     2820                    j = surf_usm_h%j(m)
     2821                    k = surf_usm_h%k(m)
     2822                    temp_pf(k,j,i) = surf_usm_h%t_surf_av(m)
     2823                 ENDDO
     2824                 DO  l = 0, 3
     2825                    DO  m = 1, surf_usm_v(l)%ns
     2826                       i = surf_usm_v(l)%i(m)
     2827                       j = surf_usm_v(l)%j(m)
     2828                       k = surf_usm_v(l)%k(m)
     2829                       temp_pf(k,j,i) = surf_usm_v(l)%t_surf_av(m)
     2830                    ENDDO
     2831                 ENDDO
     2832              ENDIF
    23902833             
    23912834          CASE ( 'usm_t_wall' )
    23922835!--           wall temperature for  iwl layer of walls and land
    2393               DO isurf = dirstart(ids), dirend(ids)
    2394                  IF ( surfl(id,isurf) == ids )  THEN
    2395                    IF ( av == 0 )  THEN
    2396                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = t_wall(iwl,isurf)
    2397                    ELSE
    2398                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = t_wall_av(iwl,isurf)
    2399                    ENDIF
    2400                  ENDIF
    2401               ENDDO
     2836              IF ( av == 0 )  THEN
     2837                 DO  m = 1, surf_usm_h%ns
     2838                    i = surf_usm_h%i(m)
     2839                    j = surf_usm_h%j(m)
     2840                    k = surf_usm_h%k(m)
     2841                    temp_pf(k,j,i) = t_wall_h(iwl,m)
     2842                 ENDDO
     2843                 DO  l = 0, 3
     2844                    DO  m = 1, surf_usm_v(l)%ns
     2845                       i = surf_usm_v(l)%i(m)
     2846                       j = surf_usm_v(l)%j(m)
     2847                       k = surf_usm_v(l)%k(m)
     2848                       temp_pf(k,j,i) = t_wall_v(l)%t(iwl,m)
     2849                    ENDDO
     2850                 ENDDO
     2851              ELSE
     2852                 DO  m = 1, surf_usm_h%ns
     2853                    i = surf_usm_h%i(m)
     2854                    j = surf_usm_h%j(m)
     2855                    k = surf_usm_h%k(m)
     2856                    temp_pf(k,j,i) = surf_usm_h%t_wall_av(iwl,m)
     2857                 ENDDO
     2858                 DO  l = 0, 3
     2859                    DO  m = 1, surf_usm_v(l)%ns
     2860                       i = surf_usm_v(l)%i(m)
     2861                       j = surf_usm_v(l)%j(m)
     2862                       k = surf_usm_v(l)%k(m)
     2863                       temp_pf(k,j,i) = surf_usm_v(l)%t_wall_av(iwl,m)
     2864                    ENDDO
     2865                 ENDDO
     2866              ENDIF
    24022867             
    24032868          CASE DEFAULT
     
    24082873!--     fill out array local_pf which is subsequently treated by data_output_3d
    24092874        CALL exchange_horiz( temp_pf, nbgp )
     2875!
     2876!--  To Do: why reversed loop order
    24102877        DO j = nysg,nyng
    24112878            DO i = nxlg,nxrg
     
    25893056        IMPLICIT NONE
    25903057
    2591         INTEGER(iwp) ::  k, l            !< running indices
     3058        INTEGER(iwp) ::  k, l, m            !< running indices
    25923059       
    25933060        CALL location_message( '    initialization of wall surface model', .TRUE. )
     
    25963063!--     Temperature is defined at the center of the wall layers,
    25973064!--     whereas gradients/fluxes are defined at the edges (_stag)
    2598         DO l = nzb_wall, nzt_wall
    2599            zwn(l) = zwn_default(l)
     3065        DO k = nzb_wall, nzt_wall
     3066           zwn(k) = zwn_default(k)
    26003067        ENDDO
    2601        
    2602 !--     apply for all particular wall grids
    2603         DO l = startenergy, endenergy
    2604            zw(:,l) = zwn(:) * thickness_wall(l)
    2605            dz_wall(nzb_wall,l) = zw(nzb_wall,l)
     3068!       
     3069!--     apply for all particular surface grids. First for horizontal surfaces
     3070        DO  m = 1, surf_usm_h%ns
     3071           surf_usm_h%zw(:,m)             = zwn(:) *                           &
     3072                                            surf_usm_h%thickness_wall(m)
     3073           surf_usm_h%dz_wall(nzb_wall,m) = surf_usm_h%zw(nzb_wall,m)
    26063074           DO k = nzb_wall+1, nzt_wall
    2607                dz_wall(k,l) = zw(k,l) - zw(k-1,l)
     3075               surf_usm_h%dz_wall(k,m) = surf_usm_h%zw(k,m) -                  &
     3076                                         surf_usm_h%zw(k-1,m)
    26083077           ENDDO
    26093078           
    2610            dz_wall(nzt_wall+1,l) = dz_wall(nzt_wall,l)
     3079           surf_usm_h%dz_wall(nzt_wall+1,m) = surf_usm_h%dz_wall(nzt_wall,m)
    26113080
    26123081           DO k = nzb_wall, nzt_wall-1
    2613                dz_wall_stag(k,l) = 0.5 * (dz_wall(k+1,l) + dz_wall(k,l))
     3082               surf_usm_h%dz_wall_stag(k,m) = 0.5 * (                          &
     3083                           surf_usm_h%dz_wall(k+1,m) + surf_usm_h%dz_wall(k,m) )
    26143084           ENDDO
    2615            dz_wall_stag(nzt_wall,l) = dz_wall(nzt_wall,l)
     3085           surf_usm_h%dz_wall_stag(nzt_wall,m) = surf_usm_h%dz_wall(nzt_wall,m)
    26163086        ENDDO
    2617        
    2618         ddz_wall      = 1.0_wp / dz_wall
    2619         ddz_wall_stag = 1.0_wp / dz_wall_stag
     3087        surf_usm_h%ddz_wall      = 1.0_wp / surf_usm_h%dz_wall
     3088        surf_usm_h%ddz_wall_stag = 1.0_wp / surf_usm_h%dz_wall_stag
     3089!       
     3090!--     For vertical surfaces
     3091        DO  l = 0, 3
     3092           DO  m = 1, surf_usm_v(l)%ns
     3093              surf_usm_v(l)%zw(:,m)             = zwn(:) *                     &
     3094                                                  surf_usm_v(l)%thickness_wall(m)
     3095              surf_usm_v(l)%dz_wall(nzb_wall,m) = surf_usm_v(l)%zw(nzb_wall,m)
     3096              DO k = nzb_wall+1, nzt_wall
     3097                  surf_usm_v(l)%dz_wall(k,m) = surf_usm_v(l)%zw(k,m) -         &
     3098                                               surf_usm_v(l)%zw(k-1,m)
     3099              ENDDO
     3100           
     3101              surf_usm_v(l)%dz_wall(nzt_wall+1,m) =                            &
     3102                                              surf_usm_v(l)%dz_wall(nzt_wall,m)
     3103
     3104              DO k = nzb_wall, nzt_wall-1
     3105                  surf_usm_v(l)%dz_wall_stag(k,m) = 0.5 * (                    &
     3106                                                surf_usm_v(l)%dz_wall(k+1,m) + &
     3107                                                surf_usm_v(l)%dz_wall(k,m) )
     3108              ENDDO
     3109              surf_usm_v(l)%dz_wall_stag(nzt_wall,m) =                         &
     3110                                              surf_usm_v(l)%dz_wall(nzt_wall,m)
     3111           ENDDO
     3112           surf_usm_v(l)%ddz_wall      = 1.0_wp / surf_usm_v(l)%dz_wall
     3113           surf_usm_v(l)%ddz_wall_stag = 1.0_wp / surf_usm_v(l)%dz_wall_stag
     3114        ENDDO     
     3115
    26203116       
    26213117        CALL location_message( '    wall structures filed out', .TRUE. )
     
    26353131        IMPLICIT NONE
    26363132
    2637         INTEGER(iwp) ::  i, j, k, l            !< running indices
     3133        INTEGER(iwp) ::  i, j, k, l, m            !< running indices
    26383134        REAL(wp)     ::  c, d, tin, exn
    26393135       
     
    26993195        ELSE
    27003196       
    2701 !--         Calculate initial surface temperature
     3197!--         Calculate initial surface temperature from pt of adjacent gridbox
    27023198            exn = ( surface_pressure / 1000.0_wp )**0.286_wp
    27033199
    2704             DO l = startenergy, endenergy
    2705                 k = surfl(iz,l)
    2706                 j = surfl(iy,l)
    2707                 i = surfl(ix,l)
    2708 
    2709 !--              Initial surface temperature set from pt of adjacent gridbox
    2710                 t_surf(l) = pt(k,j,i) * exn
     3200!
     3201!--         At horizontal surfaces. Please note, t_surf_h is defined on a
     3202!--         different data type, but with the same dimension.
     3203            DO  m = 1, surf_usm_h%ns
     3204               i = surf_usm_h%i(m)           
     3205               j = surf_usm_h%j(m)
     3206               k = surf_usm_h%k(m)
     3207
     3208               t_surf_h(m) = pt(k,j,i) * exn
    27113209            ENDDO
     3210!
     3211!--         At vertical surfaces.
     3212            DO  l = 0, 3
     3213               DO  m = 1, surf_usm_v(l)%ns
     3214                  i = surf_usm_v(l)%i(m)           
     3215                  j = surf_usm_v(l)%j(m)
     3216                  k = surf_usm_v(l)%k(m)
     3217
     3218                  t_surf_v(l)%t(m) = pt(k,j,i) * exn
     3219               ENDDO
     3220            ENDDO
     3221
    27123222     
    27133223!--         initial values for t_wall
    27143224!--         outer value is set to surface temperature
    27153225!--         inner value is set to wall_inner_temperature
    2716 !--         and profile is logaritmic (linear in nz)
    2717             DO l = startenergy, endenergy
    2718                 IF ( isroof_surf(l) )  THEN
    2719                     tin = roof_inner_temperature
    2720                 ELSE IF ( surf(id,l)==iroof )  THEN
    2721                     tin = soil_inner_temperature
    2722                 ELSE
    2723                     tin = wall_inner_temperature
    2724                 ENDIF
    2725                 DO k = nzb_wall, nzt_wall+1
    2726                     c = REAL(k-nzb_wall,wp)/REAL(nzt_wall+1-nzb_wall,wp)
    2727                     t_wall(k,:) = (1.0_wp-c)*t_surf(:) + c*tin
    2728                 ENDDO
     3226!--         and profile is logaritmic (linear in nz).
     3227!--         Horizontal surfaces
     3228            DO  m = 1, surf_usm_h%ns
     3229!
     3230!--            Roof
     3231               IF ( surf_usm_h%isroof_surf(m) )  THEN
     3232                   tin = roof_inner_temperature
     3233!
     3234!--            Normal land surface
     3235               ELSE
     3236                   tin = soil_inner_temperature
     3237               ENDIF
     3238
     3239               DO k = nzb_wall, nzt_wall+1
     3240                   c = REAL( k - nzb_wall, wp ) /                              &
     3241                       REAL( nzt_wall + 1 - nzb_wall , wp )
     3242
     3243                   t_wall_h(k,m) = ( 1.0_wp - c ) * t_surf_h(m) + c * tin
     3244               ENDDO
    27293245            ENDDO
     3246!
     3247!--         Vertical surfaces
     3248            DO  l = 0, 3
     3249               DO  m = 1, surf_usm_v(l)%ns
     3250!
     3251!--               Inner wall
     3252                  tin = wall_inner_temperature
     3253
     3254                  DO k = nzb_wall, nzt_wall+1
     3255                     c = REAL( k - nzb_wall, wp ) /                            &
     3256                         REAL( nzt_wall + 1 - nzb_wall , wp )
     3257
     3258                     t_wall_v(l)%t(k,m) = ( 1.0_wp - c ) * t_surf_v(l)%t(m) +  &
     3259                                          c * tin
     3260                  ENDDO
     3261               ENDDO
     3262            ENDDO
     3263
    27303264        ENDIF
    27313265       
    27323266!--   
    2733 !--        Possibly DO user-defined actions (e.g. define heterogeneous wall surface)
     3267!--     Possibly DO user-defined actions (e.g. define heterogeneous wall surface)
    27343268        CALL user_init_urban_surface
    27353269
    27363270!--     initialize prognostic values for the first timestep
    2737         t_surf_p = t_surf
    2738         t_wall_p = t_wall
     3271        t_surf_h_p = t_surf_h
     3272        t_surf_v_p = t_surf_v
     3273
     3274        t_wall_h_p = t_wall_h
     3275        t_wall_v_p = t_wall_v
    27393276       
    27403277!--     Adjust radiative fluxes for urban surface at model start
     
    27593296        IMPLICIT NONE
    27603297
    2761         INTEGER(iwp) ::  i,j,k,l,kw                      !< running indices
     3298        INTEGER(iwp) ::  i,j,k,l,kw, m                      !< running indices
    27623299
    27633300        REAL(wp), DIMENSION(nzb_wall:nzt_wall) :: wtend  !< tendency
    27643301
    2765                                                
    2766         DO l = startenergy, endenergy
    2767 !--         calculate frequently used parameters
    2768             k = surfl(iz,l)
    2769             j = surfl(iy,l)
    2770             i = surfl(ix,l)
    2771 
    2772             !
    2773 !--         prognostic equation for ground/wall/roof temperature t_wall
    2774             wtend(:) = 0.0_wp
    2775             wtend(nzb_wall) = (1.0_wp/rho_c_wall(nzb_wall,l)) *                     &
    2776                        ( lambda_h(nzb_wall,l) * ( t_wall(nzb_wall+1,l)              &
    2777                          - t_wall(nzb_wall,l) ) * ddz_wall(nzb_wall+1,l)            &
    2778                          + wghf_eb(l) ) * ddz_wall_stag(nzb_wall,l)
     3302!
     3303!--     For horizontal surfaces                                   
     3304        DO  m = 1, surf_usm_h%ns
     3305!
     3306!--        Obtain indices
     3307           i = surf_usm_h%i(m)           
     3308           j = surf_usm_h%j(m)
     3309           k = surf_usm_h%k(m)
     3310!
     3311!--        prognostic equation for ground/roof temperature t_wall_h
     3312           wtend(:) = 0.0_wp
     3313           wtend(nzb_wall) = (1.0_wp / surf_usm_h%rho_c_wall(nzb_wall,m)) *    &
     3314                                      ( surf_usm_h%lambda_h(nzb_wall,m) *      &
     3315                                        ( t_wall_h(nzb_wall+1,m)               &
     3316                                        - t_wall_h(nzb_wall,m) ) *             &
     3317                                        surf_usm_h%ddz_wall(nzb_wall+1,m)      &
     3318                                      + surf_usm_h%wghf_eb(m) ) *              &
     3319                                        surf_usm_h%ddz_wall_stag(nzb_wall,m)
    27793320           
    2780             DO  kw = nzb_wall+1, nzt_wall
    2781                 wtend(kw) = (1.0_wp/rho_c_wall(kw,l))                               &
    2782                               * (   lambda_h(kw,l)                                  &
    2783                                  * ( t_wall(kw+1,l) - t_wall(kw,l) )                &
    2784                                  * ddz_wall(kw+1,l)                                 &
    2785                               - lambda_h(kw-1,l)                                    &
    2786                                  * ( t_wall(kw,l) - t_wall(kw-1,l) )                &
    2787                                  * ddz_wall(kw,l)                                   &
    2788                               ) * ddz_wall_stag(kw,l)
     3321           DO  kw = nzb_wall+1, nzt_wall
     3322               wtend(kw) = (1.0_wp / surf_usm_h%rho_c_wall(kw,m))              &
     3323                              * (   surf_usm_h%lambda_h(kw,m)                  &
     3324                                 * ( t_wall_h(kw+1,m) - t_wall_h(kw,m) )       &
     3325                                 * surf_usm_h%ddz_wall(kw+1,m)                 &
     3326                              - surf_usm_h%lambda_h(kw-1,m)                    &
     3327                                 * ( t_wall_h(kw,m) - t_wall_h(kw-1,m) )       &
     3328                                 * surf_usm_h%ddz_wall(kw,m)                   &
     3329                              ) * surf_usm_h%ddz_wall_stag(kw,m)
    27893330            ENDDO
    27903331
    2791             t_wall_p(nzb_wall:nzt_wall,l) = t_wall(nzb_wall:nzt_wall,l)             &
    2792                                              + dt_3d * ( tsc(2)                     &
    2793                                              * wtend(nzb_wall:nzt_wall) + tsc(3)    &
    2794                                              * tt_wall_m(nzb_wall:nzt_wall,l) )   
     3332           t_wall_h_p(nzb_wall:nzt_wall,m) = t_wall_h(nzb_wall:nzt_wall,m)     &
     3333                                 + dt_3d * ( tsc(2)                            &
     3334                                 * wtend(nzb_wall:nzt_wall) + tsc(3)           &
     3335                                 * surf_usm_h%tt_wall_m(nzb_wall:nzt_wall,m) )   
    27953336           
    2796             !
    2797 !--         calculate t_wall tendencies for the next Runge-Kutta step
    2798             IF ( timestep_scheme(1:5) == 'runge' )  THEN
    2799                 IF ( intermediate_timestep_count == 1 )  THEN
     3337!
     3338!--        calculate t_wall tendencies for the next Runge-Kutta step
     3339           IF ( timestep_scheme(1:5) == 'runge' )  THEN
     3340               IF ( intermediate_timestep_count == 1 )  THEN
     3341                  DO  kw = nzb_wall, nzt_wall
     3342                     surf_usm_h%tt_wall_m(kw,m) = wtend(kw)
     3343                  ENDDO
     3344               ELSEIF ( intermediate_timestep_count <                          &
     3345                        intermediate_timestep_count_max )  THEN
    28003346                   DO  kw = nzb_wall, nzt_wall
    2801                       tt_wall_m(kw,l) = wtend(kw)
     3347                      surf_usm_h%tt_wall_m(kw,m) = -9.5625_wp * wtend(kw) +    &
     3348                                         5.3125_wp * surf_usm_h%tt_wall_m(kw,m)
    28023349                   ENDDO
    2803                 ELSEIF ( intermediate_timestep_count <                              &
    2804                          intermediate_timestep_count_max )  THEN
    2805                     DO  kw = nzb_wall, nzt_wall
    2806                         tt_wall_m(kw,l) = -9.5625_wp * wtend(kw) + 5.3125_wp        &
    2807                                          * tt_wall_m(kw,l)
    2808                     ENDDO
    2809                 ENDIF
    2810             ENDIF
     3350               ENDIF
     3351           ENDIF
     3352        ENDDO
     3353!
     3354!--     For vertical surfaces     
     3355        DO  l = 0, 3                             
     3356           DO  m = 1, surf_usm_v(l)%ns
     3357!
     3358!--           Obtain indices
     3359              i = surf_usm_v(l)%i(m)           
     3360              j = surf_usm_v(l)%j(m)
     3361              k = surf_usm_v(l)%k(m)
     3362!
     3363!--           prognostic equation for wall temperature t_wall_v
     3364              wtend(:) = 0.0_wp
     3365              wtend(nzb_wall) = (1.0_wp / surf_usm_v(l)%rho_c_wall(nzb_wall,m)) * &
     3366                                      ( surf_usm_v(l)%lambda_h(nzb_wall,m) *   &
     3367                                        ( t_wall_v(l)%t(nzb_wall+1,m)          &
     3368                                        - t_wall_v(l)%t(nzb_wall,m) ) *        &
     3369                                        surf_usm_v(l)%ddz_wall(nzb_wall+1,m)   &
     3370                                      + surf_usm_v(l)%wghf_eb(m) ) *           &
     3371                                        surf_usm_v(l)%ddz_wall_stag(nzb_wall,m)
     3372           
     3373              DO  kw = nzb_wall+1, nzt_wall
     3374                  wtend(kw) = (1.0_wp / surf_usm_v(l)%rho_c_wall(kw,m))        &
     3375                           * (   surf_usm_v(l)%lambda_h(kw,m)                  &
     3376                              * ( t_wall_v(l)%t(kw+1,m) - t_wall_v(l)%t(kw,m) )&
     3377                              * surf_usm_v(l)%ddz_wall(kw+1,m)                 &
     3378                           - surf_usm_v(l)%lambda_h(kw-1,m)                    &
     3379                              * ( t_wall_v(l)%t(kw,m) - t_wall_v(l)%t(kw-1,m) )&
     3380                              * surf_usm_v(l)%ddz_wall(kw,m)                   &
     3381                              ) * surf_usm_v(l)%ddz_wall_stag(kw,m)
     3382               ENDDO
     3383
     3384              t_wall_v_p(l)%t(nzb_wall:nzt_wall,m) =                           &
     3385                                   t_wall_v(l)%t(nzb_wall:nzt_wall,m)          &
     3386                                 + dt_3d * ( tsc(2)                            &
     3387                                 * wtend(nzb_wall:nzt_wall) + tsc(3)           &
     3388                                 * surf_usm_v(l)%tt_wall_m(nzb_wall:nzt_wall,m) )   
     3389           
     3390!
     3391!--           calculate t_wall tendencies for the next Runge-Kutta step
     3392              IF ( timestep_scheme(1:5) == 'runge' )  THEN
     3393                  IF ( intermediate_timestep_count == 1 )  THEN
     3394                     DO  kw = nzb_wall, nzt_wall
     3395                        surf_usm_v(l)%tt_wall_m(kw,m) = wtend(kw)
     3396                     ENDDO
     3397                  ELSEIF ( intermediate_timestep_count <                       &
     3398                           intermediate_timestep_count_max )  THEN
     3399                      DO  kw = nzb_wall, nzt_wall
     3400                         surf_usm_v(l)%tt_wall_m(kw,m) =                       &
     3401                                     - 9.5625_wp * wtend(kw) +                 &
     3402                                       5.3125_wp * surf_usm_v(l)%tt_wall_m(kw,m)
     3403                      ENDDO
     3404                  ENDIF
     3405              ENDIF
     3406           ENDDO
    28113407        ENDDO
    28123408
     
    28573453!--    Read user-defined namelist
    28583454       READ ( 11, urban_surface_par )
    2859 
    28603455!
    28613456!--    Set flag that indicates that the land surface model is switched on
    28623457       urban_surface = .TRUE.
    2863        
    28643458
    28653459 10    CONTINUE
     
    28793473        IMPLICIT NONE
    28803474       
    2881         INTEGER(iwp)               :: i, j, k, kk, is, js, d, ku, refstep
     3475        INTEGER(iwp)               :: i, j, k, kk, is, js, d, ku, refstep, m, mm, l, ll
    28823476        INTEGER(iwp)               :: nzubl, nzutl, isurf, isurfsrc, isurf1, isvf, icsf, ipcgb
    28833477        INTEGER(iwp), DIMENSION(4) :: bdycross
     
    28923486        REAL(wp)                   :: pc_box_area, pc_abs_frac, pc_abs_eff
    28933487        INTEGER(iwp)               :: pc_box_dimshift !< transform for best accuracy
     3488        INTEGER(iwp), DIMENSION(0:3) :: reorder = (/ 1, 0, 3, 2 /)
    28943489       
    28953490       
     
    29493544!--     First pass: direct + diffuse irradiance
    29503545!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    2951         surfinswdir   = 0._wp
    2952         surfinswdif   = 0._wp
    2953         surfinlwdif   = 0._wp
    2954         surfins   = 0._wp
    2955         surfinl   = 0._wp
    2956         surfoutsl    = 0._wp
    2957         surfoutll    = 0._wp
     3546        surfinswdir   = 0._wp !nsurfl
     3547        surfinswdif   = 0._wp !nsurfl
     3548        surfinlwdif   = 0._wp !nsurfl
     3549        surfins       = 0._wp !nsurfl
     3550        surfinl       = 0._wp !nsurfl
     3551        surfoutsl(:)  = 0.0_wp !start-end
     3552        surfoutll(:)  = 0.0_wp !start-end
    29583553       
    29593554!--     Set up thermal radiation from surfaces
    29603555!--     emiss_surf is defined only for surfaces for which energy balance is calculated
    2961         surfoutll(startenergy:endenergy) = emiss_surf(startenergy:endenergy) * sigma_sb   &
    2962                                            * t_surf(startenergy:endenergy)**4
     3556!--     Workaround: reorder surface data type back on 1D array including all surfaces,
     3557!--     which implies to reorder horizontal and vertical surfaces
     3558!
     3559!--     Horizontal walls
     3560        mm = 1
     3561        DO  i = nxl, nxr
     3562           DO  j = nys, nyn
     3563
     3564              DO  m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
     3565                 surfoutll(mm) = surf_usm_h%emiss_surf(m) * sigma_sb   &
     3566                                     * t_surf_h(m)**4
     3567                 albedo_surf(mm) = surf_usm_h%albedo_surf(m)
     3568                 emiss_surf(mm)  = surf_usm_h%emiss_surf(m)
     3569                 mm = mm + 1
     3570              ENDDO
     3571           ENDDO
     3572        ENDDO
     3573!
     3574!--     Vertical walls
     3575        DO  i = nxl, nxr
     3576           DO  j = nys, nyn
     3577              DO  ll = 0, 3
     3578                 l = reorder(ll)
     3579                 DO  m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i)
     3580                    surfoutll(mm) = surf_usm_v(l)%emiss_surf(m) * sigma_sb   &
     3581                                     * t_surf_v(l)%t(m)**4
     3582                    albedo_surf(mm) = surf_usm_v(l)%albedo_surf(m)
     3583                    emiss_surf(mm) = surf_usm_v(l)%emiss_surf(m)
     3584                    mm = mm + 1
     3585                 ENDDO
     3586              ENDDO
     3587           ENDDO
     3588        ENDDO
    29633589       
    29643590#if defined( __parallel )
    29653591!--     might be optimized and gather only values relevant for current processor
     3592       
    29663593        CALL MPI_AllGatherv(surfoutll, nenergy, MPI_REAL, &
    2967                             surfoutl, nsurfs, surfstart, MPI_REAL, comm2d, ierr)
     3594                            surfoutl, nsurfs, surfstart, MPI_REAL, comm2d, ierr) !nsurf global
    29683595#else
    2969         surfoutl(:) = surfoutll(:)
     3596        surfoutl(:) = surfoutll(:) !nsurf global
    29703597#endif
    29713598       
     
    30053632                surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * surfoutl(isurfsrc)
    30063633            ENDIF
    3007            
    3008             IF ( zenith(0) > 0  .AND.  all( surf(:, isurfsrc) == bdycross ) )  THEN
     3634
     3635            IF ( zenith(0) > 0  .AND.  all( surf(1:4,isurfsrc) == bdycross ) )  THEN
    30093636!--             found svf between model boundary and the face => face isn't shaded
    3010                 surfinswdir(isurf) = rad_sw_in_dir(j, i) &
     3637                surfinswdir(isurf) = rad_sw_in_dir(j,i) &
    30113638                    * costheta(surfl(id, isurf)) * svf(2,isvf) / zenith(0)
    30123639
     
    30503677                ENDIF
    30513678
    3052                 IF ( zenith(0) > 0  .AND.  all( surf(:, isurfsrc) == bdycross ) )  THEN
     3679                IF ( zenith(0) > 0  .AND.  all( surf(1:4,isurfsrc) == bdycross ) )  THEN
    30533680!--                 found svf between model boundary and the pcgb => pcgb isn't shaded
    30543681                    pc_abs_frac = 1._wp - exp(pc_abs_eff * lad_s(k,j,i))
     
    30583685            ENDDO
    30593686        ENDIF
     3687
    30603688        surfins(startenergy:endenergy) = surfinswdir(startenergy:endenergy) + surfinswdif(startenergy:endenergy)
    30613689        surfinl(startenergy:endenergy) = surfinl(startenergy:endenergy) + surfinlwdif(startenergy:endenergy)
     
    30643692        surfoutsw(:) = 0.0_wp
    30653693        surfoutlw(:) = surfoutll(:)
    3066         surfhf(startenergy:endenergy) = surfinsw(startenergy:endenergy) + surfinlw(startenergy:endenergy) &
    3067                                       - surfoutsw(startenergy:endenergy) - surfoutlw(startenergy:endenergy)
     3694!         surfhf(startenergy:endenergy) = surfinsw(startenergy:endenergy) + surfinlw(startenergy:endenergy) &
     3695!                                       - surfoutsw(startenergy:endenergy) - surfoutlw(startenergy:endenergy)
    30683696       
    30693697!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     
    31173745            surfoutsw(startenergy:endenergy) = surfoutsw(startenergy:endenergy) + surfoutsl(startenergy:endenergy)
    31183746            surfoutlw(startenergy:endenergy) = surfoutlw(startenergy:endenergy) + surfoutll(startenergy:endenergy)
    3119             surfhf(startenergy:endenergy) = surfinsw(startenergy:endenergy) + surfinlw(startenergy:endenergy) &
    3120                                           - surfoutsw(startenergy:endenergy) - surfoutlw(startenergy:endenergy)
     3747!             surfhf(startenergy:endenergy) = surfinsw(startenergy:endenergy) + surfinlw(startenergy:endenergy) &
     3748!                                           - surfoutsw(startenergy:endenergy) - surfoutlw(startenergy:endenergy)
    31213749       
    31223750        ENDDO
     
    31293757                i = pcbl(ix, ipcgb)
    31303758                k = pcbl(iz, ipcgb)
    3131                 kk = k - nzb_s_inner(j,i)  !- lad arrays are defined flat
     3759                kk = k - MAXLOC(                                               &
     3760                                MERGE( 1, 0,                                   &
     3761                                       BTEST( wall_flags_0(:,j,i), 12 )        &
     3762                                     ), DIM = 1                                &
     3763                               ) - 1 ! kk = k - nzb_s_inner(j,i)  !- lad arrays are defined flat
    31323764                pc_heating_rate(kk, j, i) = (pcbinsw(ipcgb) + pcbinlw(ipcgb)) &
    31333765                    * pchf_prep(k) * pt(k, j, i) !-- = dT/dt
    31343766            ENDDO
    31353767        ENDIF
     3768!
     3769!--     Transfer radiation arrays required for energy balance to the respective data types
     3770        DO  i = startenergy, endenergy
     3771           m  = surfl(5,i)         
     3772!
     3773!--        upward-facing
     3774           IF ( surfl(1,i) == 0 )  THEN
     3775              surf_usm_h%rad_in_sw(m)  = surfinsw(i)
     3776              surf_usm_h%rad_out_sw(m) = surfoutsw(i)
     3777              surf_usm_h%rad_in_lw(m)  = surfinlw(i)
     3778              surf_usm_h%rad_out_lw(m) = surfoutlw(i)
     3779!
     3780!--        southward-facding
     3781           ELSEIF ( surfl(1,i) == 1 )  THEN
     3782              surf_usm_v(1)%rad_in_sw(m)  = surfinsw(i)
     3783              surf_usm_v(1)%rad_out_sw(m) = surfoutsw(i)
     3784              surf_usm_v(1)%rad_in_lw(m)  = surfinlw(i)
     3785              surf_usm_v(1)%rad_out_lw(m) = surfoutlw(i)
     3786!
     3787!--        northward-facding
     3788           ELSEIF ( surfl(1,i) == 2 )  THEN
     3789              surf_usm_v(0)%rad_in_sw(m)  = surfinsw(i)
     3790              surf_usm_v(0)%rad_out_sw(m) = surfoutsw(i)
     3791              surf_usm_v(0)%rad_in_lw(m)  = surfinlw(i)
     3792              surf_usm_v(0)%rad_out_lw(m) = surfoutlw(i)
     3793!
     3794!--        westward-facding
     3795           ELSEIF ( surfl(1,i) == 3 )  THEN
     3796              surf_usm_v(3)%rad_in_sw(m)  = surfinsw(i)
     3797              surf_usm_v(3)%rad_out_sw(m) = surfoutsw(i)
     3798              surf_usm_v(3)%rad_in_lw(m)  = surfinlw(i)
     3799              surf_usm_v(3)%rad_out_lw(m) = surfoutlw(i)
     3800!
     3801!--        eastward-facing
     3802           ELSEIF ( surfl(1,i) == 4 )  THEN
     3803              surf_usm_v(2)%rad_in_sw(m)  = surfinsw(i)
     3804              surf_usm_v(2)%rad_out_sw(m) = surfoutsw(i)
     3805              surf_usm_v(2)%rad_in_lw(m)  = surfinlw(i)
     3806              surf_usm_v(2)%rad_out_lw(m) = surfoutlw(i)
     3807           ENDIF
     3808
     3809        ENDDO
     3810
     3811
     3812        DO  m = 1, surf_usm_h%ns
     3813           surf_usm_h%surfhf(m) = surf_usm_h%rad_in_sw(m)  +                   &
     3814                                  surf_usm_h%rad_in_lw(m)  -                   &
     3815                                  surf_usm_h%rad_out_sw(m) -                   &
     3816                                  surf_usm_h%rad_out_lw(m)
     3817        ENDDO
     3818
     3819        DO  l = 0, 3
     3820           DO  m = 1, surf_usm_v(l)%ns
     3821              surf_usm_v(l)%surfhf(m) = surf_usm_v(l)%rad_in_sw(m)  +          &
     3822                                        surf_usm_v(l)%rad_in_lw(m)  -          &
     3823                                        surf_usm_v(l)%rad_out_sw(m) -          &
     3824                                        surf_usm_v(l)%rad_out_lw(m)
     3825           ENDDO
     3826        ENDDO
    31363827
    31373828!--     return surface radiation to horizontal surfaces
     
    32073898        REAL(wp), PARAMETER                    :: grow_factor = 1.5_wp !< factor of expansion of grow arrays
    32083899
    3209 
     3900!
    32103901!--     Maximum number of gridboxes visited equals to maximum number of boundaries crossed in each dimension plus one. That's also
    32113902!--     the maximum number of plant canopy boxes written. We grow the acsf array accordingly using exponential factor.
     
    33514042        visible = .TRUE.
    33524043
    3353        
    33544044    END SUBROUTINE usm_raytrace
    33554045   
     
    34704160                SELECT CASE ( TRIM( variable_chr ) )
    34714161               
    3472                    CASE ( 't_surf' )
     4162                   CASE ( 't_surf_h' )
    34734163#if defined( __nopointer )                   
    3474                       IF ( .NOT.  ALLOCATED( t_surf ) )                         &
    3475                          ALLOCATE( t_surf(startenergy:endenergy) )
    3476                       READ ( 13 )  t_surf
     4164                      IF ( .NOT.  ALLOCATED( t_surf_h ) )                      &
     4165                         ALLOCATE( t_surf_h(1:surf_usm_h%ns) )
     4166                      READ ( 13 )  t_surf_h
    34774167#else                     
    3478                       IF ( .NOT.  ALLOCATED( t_surf_1 ) )                         &
    3479                          ALLOCATE( t_surf_1(startenergy:endenergy) )
    3480                       READ ( 13 )  t_surf_1
     4168                      IF ( .NOT.  ALLOCATED( t_surf_h_1 ) )                    &
     4169                         ALLOCATE( t_surf_h_1(1:surf_usm_h%ns) )
     4170                      READ ( 13 )  t_surf_h_1
    34814171#endif
    3482 
    3483                    CASE ( 't_wall' )
     4172                   CASE ( 't_surf_v(0)' )
     4173#if defined( __nopointer )                   
     4174                      IF ( .NOT.  ALLOCATED( t_surf_v(0)%t ) )                 &
     4175                         ALLOCATE( t_surf_v(0)%t(1:surf_usm_v(0)%ns) )
     4176                      READ ( 13 )  t_surf_v(0)%t
     4177#else                     
     4178                      IF ( .NOT.  ALLOCATED( t_surf_v_1(0)%t ) )               &
     4179                         ALLOCATE( t_surf_v_1(0)%t(1:surf_usm_v(0)%ns) )
     4180                      READ ( 13 )  t_surf_v_1(0)%t
     4181#endif
     4182                   CASE ( 't_surf_v(1)' )
     4183#if defined( __nopointer )                   
     4184                      IF ( .NOT.  ALLOCATED( t_surf_v(1)%t ) )                 &
     4185                         ALLOCATE( t_surf_v(1)%t(1:surf_usm_v(1)%ns) )
     4186                      READ ( 13 )  t_surf_v(1)%t
     4187#else                     
     4188                      IF ( .NOT.  ALLOCATED( t_surf_v_1(1)%t ) )               &
     4189                         ALLOCATE( t_surf_v_1(1)%t(1:surf_usm_v(1)%ns) )
     4190                      READ ( 13 )  t_surf_v_1(1)%t
     4191#endif
     4192                   CASE ( 't_surf_v(2)' )
     4193#if defined( __nopointer )                   
     4194                      IF ( .NOT.  ALLOCATED( t_surf_v(2)%t ) )                 &
     4195                         ALLOCATE( t_surf_v(2)%t(1:surf_usm_v(2)%ns) )
     4196                      READ ( 13 )  t_surf_v(2)%t
     4197#else                     
     4198                      IF ( .NOT.  ALLOCATED( t_surf_v_1(2)%t ) )               &
     4199                         ALLOCATE( t_surf_v_1(2)%t(1:surf_usm_v(2)%ns) )
     4200                      READ ( 13 )  t_surf_v_1(2)%t
     4201#endif
     4202                   CASE ( 't_surf_v(3)' )
     4203#if defined( __nopointer )                   
     4204                      IF ( .NOT.  ALLOCATED( t_surf_v(3)%t ) )                 &
     4205                         ALLOCATE( t_surf_v(3)%t(1:surf_usm_v(3)%ns) )
     4206                      READ ( 13 )  t_surf_v(3)%t
     4207#else                     
     4208                      IF ( .NOT.  ALLOCATED( t_surf_v_1(3)%t ) )               &
     4209                         ALLOCATE( t_surf_v_1(3)%t(1:surf_usm_v(3)%ns) )
     4210                      READ ( 13 )  t_surf_v_1(3)%t
     4211#endif
     4212                   CASE ( 't_wall_h' )
    34844213#if defined( __nopointer )
    3485                       IF ( .NOT.  ALLOCATED( t_wall ) )                         &
    3486                          ALLOCATE( t_wall(nzb_wall:nzt_wall+1,startenergy:endenergy) )
    3487                       READ ( 13 )  t_wall
     4214                      IF ( .NOT.  ALLOCATED( t_wall_h ) )                      &
     4215                         ALLOCATE( t_wall_h(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
     4216                      READ ( 13 )  t_wall_h
    34884217#else
    3489                       IF ( .NOT.  ALLOCATED( t_wall_1 ) )                         &
    3490                          ALLOCATE( t_wall_1(nzb_wall:nzt_wall+1,startenergy:endenergy) )
    3491                       READ ( 13 )  t_wall_1
     4218                      IF ( .NOT.  ALLOCATED( t_wall_h_1 ) )                    &
     4219                         ALLOCATE( t_wall_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
     4220                      READ ( 13 )  t_wall_h_1
     4221#endif
     4222                   CASE ( 't_wall_v(0)' )
     4223#if defined( __nopointer )
     4224                      IF ( .NOT.  ALLOCATED( t_wall_v(0)%t ) )                      &
     4225                         ALLOCATE( t_wall_v(0)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(0)%ns) )
     4226                      READ ( 13 )  t_wall_v(0)%t
     4227#else
     4228                      IF ( .NOT.  ALLOCATED( t_wall_v_1(0)%t ) )                    &
     4229                         ALLOCATE( t_wall_v_1(0)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(0)%ns) )
     4230                      READ ( 13 )  t_wall_v_1(0)%t
     4231#endif
     4232                   CASE ( 't_wall_v(1)' )
     4233#if defined( __nopointer )
     4234                      IF ( .NOT.  ALLOCATED( t_wall_v(1)%t ) )                      &
     4235                         ALLOCATE( t_wall_v(1)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(1)%ns) )
     4236                      READ ( 13 )  t_wall_v(1)%t
     4237#else
     4238                      IF ( .NOT.  ALLOCATED( t_wall_v_1(0)%t ) )                    &
     4239                         ALLOCATE( t_wall_v_1(1)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(1)%ns) )
     4240                      READ ( 13 )  t_wall_v_1(1)%t
     4241#endif
     4242                   CASE ( 't_wall_v(2)' )
     4243#if defined( __nopointer )
     4244                      IF ( .NOT.  ALLOCATED( t_wall_v(2)%t ) )                      &
     4245                         ALLOCATE( t_wall_v(2)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(2)%ns) )
     4246                      READ ( 13 )  t_wall_v(2)%t
     4247#else
     4248                      IF ( .NOT.  ALLOCATED( t_wall_v_1(2)%t ) )                    &
     4249                         ALLOCATE( t_wall_v_1(2)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(2)%ns) )
     4250                      READ ( 13 )  t_wall_v_1(2)%t
     4251#endif
     4252                   CASE ( 't_wall_v(3)' )
     4253#if defined( __nopointer )
     4254                      IF ( .NOT.  ALLOCATED( t_wall_v(3)%t ) )                      &
     4255                         ALLOCATE( t_wall_v(3)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(3)%ns) )
     4256                      READ ( 13 )  t_wall_v(3)%t
     4257#else
     4258                      IF ( .NOT.  ALLOCATED( t_wall_v_1(3)%t ) )                    &
     4259                         ALLOCATE( t_wall_v_1(3)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(3)%ns) )
     4260                      READ ( 13 )  t_wall_v_1(3)%t
    34924261#endif
    34934262
     
    35934362        INTEGER(iwp), DIMENSION(0:17, nysg:nyng, nxlg:nxrg)   :: usm_par
    35944363        REAL(wp), DIMENSION(1:14, nysg:nyng, nxlg:nxrg)       :: usm_val
    3595         INTEGER(iwp)                                          :: k, l, d, iw, jw, kw, it, ip, ii, ij
     4364        INTEGER(iwp)                                          :: k, l, d, iw, jw, kw, it, ip, ii, ij, m
    35964365        INTEGER(iwp)                                          :: i, j
    35974366        INTEGER(iwp)                                          :: nz, roof, dirwe, dirsn
     
    37334502            ENDDO
    37344503        ENDDO
    3735        
    3736 !--     assign the surface types to local surface array
    3737         DO  l = startenergy, endenergy
    3738            
    3739             d = surfl(id,l)
    3740             kw = surfl(iz,l)
    3741             j = surfl(iy,l)
    3742             i = surfl(ix,l)
    3743             IF ( d == iroof )  THEN
    3744 !--             horizontal surface - land or roof
    3745                 iw = i
    3746                 jw = j
    3747                 IF ( usm_par(5,jw,iw) == 0 )  THEN
    3748                     IF ( zu(kw) >= roof_height_limit )  THEN
    3749                         isroof_surf(l) = .TRUE.
    3750                         surface_types(l) = roof_category         !< default category for root surface
    3751                     ELSE
    3752                         isroof_surf(l) = .FALSE.
    3753                         surface_types(l) = land_category         !< default category for land surface
    3754                     ENDIF
    3755                     albedo_surf(l) = -1.0_wp
    3756                     thickness_wall(l) = -1.0_wp
    3757                 ELSE
    3758                     IF ( usm_par(2,jw,iw)==0 )  THEN
    3759                         isroof_surf(l) = .FALSE.
    3760                         thickness_wall(l) = -1.0_wp
    3761                     ELSE
    3762                         isroof_surf(l) = .TRUE.
    3763                         thickness_wall(l) = usm_val(2,jw,iw)
    3764                     ENDIF
    3765                     surface_types(l) = usm_par(5,jw,iw)
    3766                     albedo_surf(l) = usm_val(1,jw,iw)
    3767                 ENDIF
    3768             ELSE
    3769                 SELECT CASE (d)
    3770                     CASE (iwest)
    3771                         iw = i
    3772                         jw = j
    3773                         ii = 6
    3774                         ij = 3
    3775                     CASE (ieast)
    3776                         iw = i-1
    3777                         jw = j
    3778                         ii = 6
    3779                         ij = 3
    3780                     CASE (isouth)
    3781                         iw = i
    3782                         jw = j
    3783                         ii = 12
    3784                         ij = 9
    3785                     CASE (inorth)
    3786                         iw = i
    3787                         jw = j-1
    3788                         ii = 12
    3789                         ij = 9
    3790                 END SELECT
    3791                
    3792                 IF ( kw <= usm_par(ii,jw,iw) )  THEN
    3793 !--                 pedestrant zone
    3794                     isroof_surf(l) = .FALSE.
    3795                     IF ( usm_par(ii+1,jw,iw) == 0 )  THEN
    3796                         surface_types(l) = pedestrant_category   !< default category for wall surface in pedestrant zone
    3797                         albedo_surf(l) = -1.0_wp
    3798                         thickness_wall(l) = -1.0_wp
    3799                     ELSE
    3800                         surface_types(l) = usm_par(ii+1,jw,iw)
    3801                         albedo_surf(l) = usm_val(ij,jw,iw)
    3802                         thickness_wall(l) = usm_val(ij+1,jw,iw)
    3803                     ENDIF
    3804                 ELSE IF ( kw <= usm_par(ii+2,jw,iw) )  THEN
    3805 !--                 wall zone
    3806                     isroof_surf(l) = .FALSE.
    3807                     IF ( usm_par(ii+3,jw,iw) == 0 )  THEN
    3808                         surface_types(l) = wall_category         !< default category for wall surface
    3809                         albedo_surf(l) = -1.0_wp
    3810                         thickness_wall(l) = -1.0_wp
    3811                     ELSE
    3812                         surface_types(l) = usm_par(ii+3,jw,iw)
    3813                         albedo_surf(l) = usm_val(ij+2,jw,iw)
    3814                         thickness_wall(l) = usm_val(ij+3,jw,iw)
    3815                     ENDIF
    3816                 ELSE IF ( kw <= usm_par(ii+4,jw,iw) )  THEN
    3817 !--                 roof zone
    3818                     isroof_surf(l) = .TRUE.
    3819                     IF ( usm_par(ii+5,jw,iw) == 0 )  THEN
    3820                         surface_types(l) = roof_category         !< default category for roof surface
    3821                         albedo_surf(l) = -1.0_wp
    3822                         thickness_wall(l) = -1.0_wp
    3823                     ELSE
    3824                         surface_types(l) = usm_par(ii+5,jw,iw)
    3825                         albedo_surf(l) = usm_val(ij+4,jw,iw)
    3826                         thickness_wall(l) = usm_val(ij+5,jw,iw)
    3827                     ENDIF
    3828                 ELSE
    3829 !--                 something wrong
    3830                     CALL message( 'usm_read_urban_surface', 'PA0505', 1, 2, 0, 6, 0 )
    3831                 ENDIF
    3832             ENDIF
    3833            
    3834 !--         find the type position
    3835             it = surface_types(l)
    3836             ip = -99999
    3837             DO k = 1, n_surface_types
    3838                 IF ( surface_type_codes(k) == it )  THEN
     4504!       
     4505!--     Assign the surface types to the respective data type.
     4506!--     First, for horizontal upward-facing surfaces.
     4507        DO  m = 1, surf_usm_h%ns
     4508           iw = surf_usm_h%i(m)
     4509           jw = surf_usm_h%j(m)
     4510           kw = surf_usm_h%k(m)
     4511
     4512           IF ( usm_par(5,jw,iw) == 0 )  THEN
     4513              IF ( zu(kw) >= roof_height_limit )  THEN
     4514                 surf_usm_h%isroof_surf(m)   = .TRUE.
     4515                 surf_usm_h%surface_types(m) = roof_category         !< default category for root surface
     4516              ELSE
     4517                 surf_usm_h%isroof_surf(m)   = .FALSE.
     4518                 surf_usm_h%surface_types(m) = land_category         !< default category for land surface
     4519              ENDIF
     4520              surf_usm_h%albedo_surf(m)    = -1.0_wp
     4521              surf_usm_h%thickness_wall(m) = -1.0_wp
     4522           ELSE
     4523              IF ( usm_par(2,jw,iw)==0 )  THEN
     4524                 surf_usm_h%isroof_surf(m)    = .FALSE.
     4525                 surf_usm_h%thickness_wall(m) = -1.0_wp
     4526              ELSE
     4527                 surf_usm_h%isroof_surf(m)    = .TRUE.
     4528                 surf_usm_h%thickness_wall(m) = usm_val(2,jw,iw)
     4529              ENDIF
     4530              surf_usm_h%surface_types(m) = usm_par(5,jw,iw)
     4531              surf_usm_h%albedo_surf(m)   = usm_val(1,jw,iw)
     4532           ENDIF
     4533!
     4534!--        Find the type position
     4535           it = surf_usm_h%surface_types(m)
     4536           ip = -99999
     4537           DO k = 1, n_surface_types
     4538              IF ( surface_type_codes(k) == it )  THEN
     4539                 ip = k
     4540                 EXIT
     4541              ENDIF
     4542           ENDDO
     4543           IF ( ip == -99999 )  THEN
     4544!--           wall category not found
     4545              WRITE (message_string, "(A,I5,A,3I5)") 'wall category ', it,     &
     4546                                     ' not found  for i,j,k=', iw,jw,kw
     4547              CALL message( 'usm_read_urban_surface', 'PA0506', 1, 2, 0, 6, 0 )
     4548           ENDIF
     4549!
     4550!--        Albedo
     4551           IF ( surf_usm_h%albedo_surf(m) < 0.0_wp )  THEN
     4552              surf_usm_h%albedo_surf(m) = surface_params(ialbedo,ip)
     4553           ENDIF
     4554!
     4555!--        emissivity of the wall
     4556           surf_usm_h%emiss_surf(m) = surface_params(iemiss,ip)
     4557!           
     4558!--        heat conductivity λS between air and wall ( W m−2 K−1 )
     4559           surf_usm_h%lambda_surf(m) = surface_params(ilambdas,ip)
     4560!           
     4561!--        roughness relative to concrete
     4562           surf_usm_h%roughness_wall(m) = surface_params(irough,ip)
     4563!           
     4564!--        Surface skin layer heat capacity (J m−2 K−1 )
     4565           surf_usm_h%c_surface(m) = surface_params(icsurf,ip)
     4566!           
     4567!--        wall material parameters:
     4568!--        thickness of the wall (m)
     4569!--        missing values are replaced by default value for category
     4570           IF ( surf_usm_h%thickness_wall(m) <= 0.001_wp )  THEN
     4571                surf_usm_h%thickness_wall(m) = surface_params(ithick,ip)
     4572           ENDIF
     4573!           
     4574!--        volumetric heat capacity rho*C of the wall ( J m−3 K−1 )
     4575           surf_usm_h%rho_c_wall(:,m) = surface_params(irhoC,ip)
     4576!           
     4577!--        thermal conductivity λH of the wall (W m−1 K−1 )
     4578           surf_usm_h%lambda_h(:,m) = surface_params(ilambdah,ip)
     4579
     4580        ENDDO
     4581!
     4582!--     For vertical surface elements ( 0 -- northward-facing, 1 -- southward-facing,
     4583!--     2 -- eastward-facing, 3 -- westward-facing )
     4584        DO  l = 0, 3
     4585           DO  m = 1, surf_usm_v(l)%ns
     4586              i  = surf_usm_v(l)%i(m)
     4587              j  = surf_usm_v(l)%j(m)
     4588              kw = surf_usm_v(l)%k(m)
     4589
     4590              IF ( l == 3 )  THEN ! westward facing
     4591                 iw = i
     4592                 jw = j
     4593                 ii = 6
     4594                 ij = 3
     4595              ELSEIF ( l == 2 )  THEN
     4596                 iw = i-1
     4597                 jw = j
     4598                 ii = 6
     4599                 ij = 3
     4600              ELSEIF ( l == 1 )  THEN
     4601                 iw = i
     4602                 jw = j
     4603                 ii = 12
     4604                 ij = 9
     4605              ELSEIF ( l == 0 )  THEN
     4606                 iw = i
     4607                 jw = j-1
     4608                 ii = 12
     4609                 ij = 9
     4610              ENDIF
     4611
     4612              IF ( kw <= usm_par(ii,jw,iw) )  THEN
     4613!--              pedestrant zone
     4614                 IF ( usm_par(ii+1,jw,iw) == 0 )  THEN
     4615                     surf_usm_v(l)%surface_types(m)  = pedestrant_category   !< default category for wall surface in pedestrant zone
     4616                     surf_usm_v(l)%albedo_surf(m)    = -1.0_wp
     4617                     surf_usm_v(l)%thickness_wall(m) = -1.0_wp
     4618                 ELSE
     4619                     surf_usm_v(l)%surface_types(m)  = usm_par(ii+1,jw,iw)
     4620                     surf_usm_v(l)%albedo_surf(m)    = usm_val(ij,jw,iw)
     4621                     surf_usm_v(l)%thickness_wall(m) = usm_val(ij+1,jw,iw)
     4622                 ENDIF
     4623              ELSE IF ( kw <= usm_par(ii+2,jw,iw) )  THEN
     4624!--              wall zone
     4625                 IF ( usm_par(ii+3,jw,iw) == 0 )  THEN
     4626                     surf_usm_v(l)%surface_types(m)  = wall_category         !< default category for wall surface
     4627                     surf_usm_v(l)%albedo_surf(m)    = -1.0_wp
     4628                     surf_usm_v(l)%thickness_wall(m) = -1.0_wp
     4629                 ELSE
     4630                     surf_usm_v(l)%surface_types(m)  = usm_par(ii+3,jw,iw)
     4631                     surf_usm_v(l)%albedo_surf(m)    = usm_val(ij+2,jw,iw)
     4632                     surf_usm_v(l)%thickness_wall(m) = usm_val(ij+3,jw,iw)
     4633                 ENDIF
     4634              ELSE IF ( kw <= usm_par(ii+4,jw,iw) )  THEN
     4635!--              roof zone
     4636                 IF ( usm_par(ii+5,jw,iw) == 0 )  THEN
     4637                     surf_usm_v(l)%surface_types(m)  = roof_category         !< default category for roof surface
     4638                     surf_usm_v(l)%albedo_surf(m)    = -1.0_wp
     4639                     surf_usm_v(l)%thickness_wall(m) = -1.0_wp
     4640                 ELSE
     4641                     surf_usm_v(l)%surface_types(m)  = usm_par(ii+5,jw,iw)
     4642                     surf_usm_v(l)%albedo_surf(m)    = usm_val(ij+4,jw,iw)
     4643                     surf_usm_v(l)%thickness_wall(m) = usm_val(ij+5,jw,iw)
     4644                 ENDIF
     4645              ELSE
     4646!--              something wrong
     4647                 CALL message( 'usm_read_urban_surface', 'PA0505', 1, 2, 0, 6, 0 )
     4648              ENDIF
     4649
     4650!
     4651!--           Find the type position
     4652              it = surf_usm_v(l)%surface_types(m)
     4653              ip = -99999
     4654              DO k = 1, n_surface_types
     4655                 IF ( surface_type_codes(k) == it )  THEN
    38394656                    ip = k
    38404657                    EXIT
    3841                 ENDIF
    3842             ENDDO
    3843             IF ( ip == -99999 )  THEN
    3844 !--             wall category not found
    3845                 WRITE (message_string, "(A,I5,A,3I5)") 'wall category ', it, ' not found  for i,j,k=', iw,jw,kw
    3846                 CALL message( 'usm_read_urban_surface', 'PA0506', 1, 2, 0, 6, 0 )
    3847             ENDIF
    3848            
    3849 !--         Fill out the parameters of the wall
    3850 !--         wall surface:
    3851            
    3852 !--         albedo
    3853             IF ( albedo_surf(l) < 0.0_wp )  THEN
    3854                 albedo_surf(l) = surface_params(ialbedo, ip)
    3855             ENDIF
    3856            
    3857 !--         emissivity of the wall
    3858             emiss_surf(l) = surface_params(iemiss, ip)
    3859            
    3860 !--         heat conductivity λS between air and wall ( W m−2 K−1 )
    3861             lambda_surf(l) = surface_params(ilambdas, ip)
    3862            
    3863 !--         roughness relative to concrete
    3864             roughness_wall(l) = surface_params(irough, ip)
    3865            
    3866 !--         Surface skin layer heat capacity (J m−2 K−1 )
    3867             c_surface(l) = surface_params(icsurf, ip)
    3868            
    3869 !--         wall material parameters:
    3870            
    3871 !--         thickness of the wall (m)
    3872 !--         missing values are replaced by default value for category
    3873             IF ( thickness_wall(l) <= 0.001_wp )  THEN
    3874                 thickness_wall(l) = surface_params(ithick, ip)
    3875             ENDIF
    3876            
    3877 !--         volumetric heat capacity rho*C of the wall ( J m−3 K−1 )
    3878             rho_c_wall(:,l) = surface_params(irhoC, ip)
    3879            
    3880 !--         thermal conductivity λH of the wall (W m−1 K−1 )
    3881             lambda_h(:,l) = surface_params(ilambdah, ip)
    3882            
     4658                 ENDIF
     4659              ENDDO
     4660              IF ( ip == -99999 )  THEN
     4661!--              wall category not found
     4662                 WRITE (message_string, "(A,I5,A,3I5)") 'wall category ', it,  &
     4663                                        ' not found  for i,j,k=', iw,jw,kw
     4664                 CALL message( 'usm_read_urban_surface', 'PA0506', 1, 2, 0, 6, 0 )
     4665              ENDIF
     4666!
     4667!--           Albedo
     4668              IF ( surf_usm_v(l)%albedo_surf(m) < 0.0_wp )  THEN
     4669                 surf_usm_v(l)%albedo_surf(m) = surface_params(ialbedo,ip)
     4670              ENDIF
     4671!
     4672!--           emissivity of the wall
     4673              surf_usm_v(l)%emiss_surf(m) = surface_params(iemiss,ip)
     4674!           
     4675!--           heat conductivity λS between air and wall ( W m−2 K−1 )
     4676              surf_usm_v(l)%lambda_surf(m) = surface_params(ilambdas,ip)
     4677!           
     4678!--           roughness relative to concrete
     4679              surf_usm_v(l)%roughness_wall(m) = surface_params(irough,ip)
     4680!           
     4681!--           Surface skin layer heat capacity (J m−2 K−1 )
     4682              surf_usm_v(l)%c_surface(m) = surface_params(icsurf,ip)
     4683!           
     4684!--           wall material parameters:
     4685!--           thickness of the wall (m)
     4686!--           missing values are replaced by default value for category
     4687              IF ( surf_usm_v(l)%thickness_wall(m) <= 0.001_wp )  THEN
     4688                   surf_usm_v(l)%thickness_wall(m) = surface_params(ithick,ip)
     4689              ENDIF
     4690!           
     4691!--           volumetric heat capacity rho*C of the wall ( J m−3 K−1 )
     4692              surf_usm_v(l)%rho_c_wall(:,m) = surface_params(irhoC,ip)
     4693!           
     4694!--           thermal conductivity λH of the wall (W m−1 K−1 )
     4695              surf_usm_v(l)%lambda_h(:,m) = surface_params(ilambdah,ip)
     4696
     4697           ENDDO
    38834698        ENDDO
    38844699
     
    39004715        IMPLICIT NONE
    39014716
    3902         INTEGER(iwp)                          :: i, j, k, l, d      !< running indices
     4717        INTEGER(iwp)                          :: i, j, k, l, d, m   !< running indices
    39034718       
    39044719        REAL(wp)                              :: pt1                !< temperature at first grid box adjacent to surface
     
    39244739       
    39254740        exn(:) = (hyp(nzub:nzut) / 100000.0_wp )**0.286_wp          !< Exner function
     4741!       
     4742!--     First, treat horizontal surface elements
     4743
     4744        DO  m = 1, surf_usm_h%ns
     4745!
     4746!--        Get indices of respective grid point
     4747           i = surf_usm_h%i(m)
     4748           j = surf_usm_h%j(m)
     4749           k = surf_usm_h%k(m)
     4750!
     4751!--        TODO - how to calculate lambda_surface for horizontal surfaces
     4752!--        (lambda_surface is set according to stratification in land surface model)
     4753!--        MS: ???
     4754           IF ( surf_usm_h%ol(m) >= 0.0_wp )  THEN
     4755              lambda_surface = surf_usm_h%lambda_surf(m)
     4756           ELSE
     4757              lambda_surface = surf_usm_h%lambda_surf(m)
     4758           ENDIF
    39264759           
    3927 !--   
    3928         DO l = startenergy, endenergy
    3929 !--         Calculate frequently used parameters
    3930             d = surfl(id,l)
    3931             k = surfl(iz,l)
    3932             j = surfl(iy,l)
    3933             i = surfl(ix,l)
    3934 
    3935 !--         TODO - how to calculate lambda_surface for horizontal surfaces
    3936 !--         (lambda_surface is set according to stratification in land surface model)
    3937             IF ( ol(j,i) >= 0.0_wp )  THEN
    3938                 lambda_surface = lambda_surf(l)
    3939             ELSE
    3940                 lambda_surface = lambda_surf(l)
    3941             ENDIF
     4760           pt1  = pt(k,j,i)
     4761!
     4762!--        calculate rho * cp coefficient at surface layer
     4763           rho_cp  = cp * hyp(k) / ( r_d * pt1 * exn(k) )
     4764!
     4765!--        Calculate aerodyamic resistance.
     4766!--        Calculation for horizontal surfaces follows LSM formulation
     4767!--        pt, us, ts are not available for the prognostic time step,
     4768!--        data from the last time step is used here.
     4769               
     4770           r_a = ( pt1 - t_surf_h(m) / exn(k) ) /                              &
     4771                 ( surf_usm_h%ts(m) * surf_usm_h%us(m) + 1.0E-10_wp )
     4772               
     4773!--        make sure that the resistance does not drop to zero
     4774           IF ( ABS(r_a) < 1.0E-10_wp )  r_a = 1.0E-10_wp
     4775               
     4776!--        the parameterization is developed originally for larger scales
     4777!--        (compare with remark in TUF-3D)
     4778!--        our first experiences show that the parameterization underestimates
     4779!--        r_a in meter resolution.
     4780!--        temporary solution - multiplication by magic constant :-(.
     4781           r_a = r_a * ra_horiz_coef
     4782               
     4783!--        factor for shf_eb
     4784           f_shf  = rho_cp / r_a
     4785       
     4786!--        add LW up so that it can be removed in prognostic equation
     4787           surf_usm_h%rad_net_l(m) = surf_usm_h%rad_in_sw(m)  -                &
     4788                                     surf_usm_h%rad_out_sw(m) +                &
     4789                                     surf_usm_h%rad_in_lw(m)  -                &
     4790                                     surf_usm_h%rad_out_lw(m)
     4791
     4792!--        numerator of the prognostic equation
     4793           coef_1 = surf_usm_h%rad_net_l(m) +                                  &
     4794                   ( 3.0_wp + 1.0_wp ) * surf_usm_h%emiss_surf(m) * sigma_sb * &
     4795                                       t_surf_h(m) ** 4 +                      & 
     4796                                       f_shf * pt1 +                           &
     4797                                       lambda_surface * t_wall_h(nzb_wall,m)
     4798
     4799!--        denominator of the prognostic equation
     4800           coef_2 = 4.0_wp * surf_usm_h%emiss_surf(m) * sigma_sb *            &
     4801                             t_surf_h(m) ** 3                                 &
     4802                           + lambda_surface + f_shf / exn(k)
     4803
     4804!--        implicit solution when the surface layer has no heat capacity,
     4805!--        otherwise use RK3 scheme.
     4806           t_surf_h_p(m) = ( coef_1 * dt_3d * tsc(2) +                        &
     4807                             surf_usm_h%c_surface(m) * t_surf_h(m) ) /        &
     4808                           ( surf_usm_h%c_surface(m) + coef_2 * dt_3d * tsc(2) )
     4809
     4810!--        add RK3 term
     4811           t_surf_h_p(m) = t_surf_h_p(m) + dt_3d * tsc(3) *                   &
     4812                           surf_usm_h%tt_surface_m(m)
    39424813           
    3943             pt1  = pt(k,j,i)
    3944 
    3945 !--         calculate rho * cp coefficient at surface layer
    3946             rho_cp  = cp * hyp(k) / ( r_d * pt1 * exn(k) )
    3947 
    3948 !--         calculate aerodyamic resistance.
    3949             IF ( d == iroof )  THEN
    3950 !--             calculation for horizontal surfaces follows LSM formulation
    3951 !--             pt, us, ts are not available for the prognostic time step,
    3952 !--             data from the last time step is used here.
     4814!--        calculate true tendency
     4815           stend = ( t_surf_h_p(m) - t_surf_h(m) - dt_3d * tsc(3) *           &
     4816                     surf_usm_h%tt_surface_m(m)) / ( dt_3d  * tsc(2) )
     4817
     4818!--        calculate t_surf tendencies for the next Runge-Kutta step
     4819           IF ( timestep_scheme(1:5) == 'runge' )  THEN
     4820              IF ( intermediate_timestep_count == 1 )  THEN
     4821                 surf_usm_h%tt_surface_m(m) = stend
     4822              ELSEIF ( intermediate_timestep_count <                          &
     4823                        intermediate_timestep_count_max )  THEN
     4824                 surf_usm_h%tt_surface_m(m) = -9.5625_wp * stend +            &
     4825                                     5.3125_wp * surf_usm_h%tt_surface_m(m)
     4826              ENDIF
     4827           ENDIF
     4828
     4829!--        in case of fast changes in the skin temperature, it is required to
     4830!--        update the radiative fluxes in order to keep the solution stable
     4831           IF ( ABS( t_surf_h_p(m) - t_surf_h(m) ) > 1.0_wp )  THEN
     4832              force_radiation_call_l = .TRUE.
     4833           ENDIF
     4834!           
     4835!--        for horizontal surfaces is pt(nzb_s_inner(j,i),j,i) = pt_surf.
     4836!--        there is no equivalent surface gridpoint for vertical surfaces.
     4837!--        pt(k,j,i) is calculated for all directions in diffusion_s
     4838!--        using surface and wall heat fluxes
     4839           pt(k-1,j,i) = t_surf_h_p(m) / exn(k)  ! not for vertical surfaces
     4840
     4841!--        calculate fluxes
     4842!--        rad_net_l is never used!           
     4843           surf_usm_h%rad_net_l(m) = surf_usm_h%rad_net_l(m) +                &
     4844                                     3.0_wp * sigma_sb *                      &
     4845                                     t_surf_h(m)**4 - 4.0_wp * sigma_sb *     &
     4846                                     t_surf_h(m)**3 * t_surf_h_p(m)
     4847           surf_usm_h%wghf_eb(m)   = lambda_surface *                         &
     4848                                      ( t_surf_h_p(m) - t_wall_h(nzb_wall,m) )
     4849!
     4850!--        ground/wall/roof surface heat flux
     4851           surf_usm_h%wshf_eb(m)   = - f_shf  * ( pt1 - t_surf_h_p(m) )
     4852!           
     4853!--        store kinematic surface heat fluxes for utilization in other processes
     4854!--        diffusion_s, surface_layer_fluxes,...
     4855           surf_usm_h%shf(m) = surf_usm_h%wshf_eb(m) / rho_cp
     4856
     4857       ENDDO
     4858!
     4859!--    Now, treat vertical surface elements
     4860       DO  l = 0, 3
     4861          DO  m = 1, surf_usm_v(l)%ns
     4862!
     4863!--          Get indices of respective grid point
     4864             i = surf_usm_v(l)%i(m)
     4865             j = surf_usm_v(l)%j(m)
     4866             k = surf_usm_v(l)%k(m)
     4867!
     4868!--          TODO - how to calculate lambda_surface for horizontal (??? do you mean verical ???) surfaces
     4869!--          (lambda_surface is set according to stratification in land surface model).
     4870!--          Please note, for vertical surfaces no ol is defined, since
     4871!--          stratification is not considered in this case.
     4872             lambda_surface = surf_usm_v(l)%lambda_surf(m)
     4873           
     4874             pt1  = pt(k,j,i)
     4875!
     4876!--          calculate rho * cp coefficient at surface layer
     4877             rho_cp  = cp * hyp(k) / ( r_d * pt1 * exn(k) )
     4878
     4879!--          Calculation of r_a for vertical surfaces
     4880!--
     4881!--          heat transfer coefficient for forced convection along vertical walls
     4882!--          follows formulation in TUF3d model (Krayenhoff & Voogt, 2006)
     4883!--           
     4884!--          H = httc (Tsfc - Tair)
     4885!--          httc = rw * (11.8 + 4.2 * Ueff) - 4.0
     4886!--           
     4887!--                rw: wall patch roughness relative to 1.0 for concrete
     4888!--                Ueff: effective wind speed
     4889!--                - 4.0 is a reduction of Rowley et al (1930) formulation based on
     4890!--                Cole and Sturrock (1977)
     4891!--           
     4892!--                Ucan: Canyon wind speed
     4893!--                wstar: convective velocity
     4894!--                Qs: surface heat flux
     4895!--                zH: height of the convective layer
     4896!--                wstar = (g/Tcan*Qs*zH)**(1./3.)
    39534897               
    3954                 r_a = (pt1 - t_surf(l)/exn(k)) / (ts(j,i) * us(j,i) + 1.0E-10_wp)
     4898!--          Effective velocity components must always
     4899!--          be defined at scalar grid point. The wall normal component is
     4900!--          obtained by simple linear interpolation. ( An alternative would
     4901!--          be an logarithmic interpolation. )
     4902             u1 = ( u(k,j,i) + u(k,j,i+1) ) * 0.5_wp
     4903             v1 = ( v(k,j,i) + v(k,j+1,i) ) * 0.5_wp
     4904             w1 = ( w(k,j,i) + w(k-1,j,i) ) * 0.5_wp
    39554905               
    3956 !--             make sure that the resistance does not drop to zero
    3957                 IF ( ABS(r_a) < 1.0E-10_wp )  r_a = 1.0E-10_wp
    3958                
    3959 !--             the parameterization is developed originally for larger scales
    3960 !--             (compare with remark in TUF-3D)
    3961 !--             our first experiences show that the parameterization underestimates
    3962 !--             r_a in meter resolution.
    3963 !--             temporary solution - multiplication by magic constant :-(.
    3964                 r_a = r_a * ra_horiz_coef
    3965                
    3966 !--             factor for shf_eb
    3967                 f_shf  = rho_cp / r_a
    3968             ELSE
    3969 !--             calculation of r_a for vertical surfaces
    3970 !--
    3971 !--             heat transfer coefficient for forced convection along vertical walls
    3972 !--             follows formulation in TUF3d model (Krayenhoff & Voogt, 2006)
    3973 !--           
    3974 !--             H = httc (Tsfc - Tair)
    3975 !--             httc = rw * (11.8 + 4.2 * Ueff) - 4.0
    3976 !--           
    3977 !--                   rw: wall patch roughness relative to 1.0 for concrete
    3978 !--                   Ueff: effective wind speed
    3979 !--                   - 4.0 is a reduction of Rowley et al (1930) formulation based on
    3980 !--                   Cole and Sturrock (1977)
    3981 !--           
    3982 !--                   Ucan: Canyon wind speed
    3983 !--                   wstar: convective velocity
    3984 !--                   Qs: surface heat flux
    3985 !--                   zH: height of the convective layer
    3986 !--                   wstar = (g/Tcan*Qs*zH)**(1./3.)
    3987                
    3988 !--             staggered grid needs to be taken into consideration
    3989                 IF ( d == inorth )  THEN
    3990                     u1 = (u(k,j,i)+u(k,j,i+1))*0.5_wp
    3991                     v1 = v(k,j+1,i)
    3992                 ELSE IF ( d == isouth )  THEN
    3993                     u1 = (u(k,j,i)+u(k,j,i+1))*0.5_wp
    3994                     v1 = v(k,j,i)
    3995                 ELSE IF ( d == ieast )  THEN
    3996                     u1 = u(k,j,i+1)
    3997                     v1 = (v(k,j,i)+v(k,j+1,i))*0.5_wp
    3998                 ELSE IF ( d == iwest )  THEN
    3999                     u1 = u(k,j,i)
    4000                     v1 = (v(k,j,i)+v(k,j+1,i))*0.5_wp
    4001                 ELSE
    4002                     STOP
    4003                 ENDIF
    4004                 w1 = (w(k,j,i)+w(k-1,j,i))*0.5_wp
    4005                
    4006                 Ueff = SQRT(u1**2 + v1**2 + w1**2)
    4007                 httc = roughness_wall(l) * (11.8 + 4.2 * Ueff) - 4.0
    4008                 f_shf  = httc
    4009             ENDIF
    4010        
    4011 !--         add LW up so that it can be removed in prognostic equation
    4012             rad_net_l(l) = surfinsw(l) - surfoutsw(l) + surfinlw(l) - surfoutlw(l)
    4013 
    4014 !--         numerator of the prognostic equation
    4015             coef_1 = rad_net_l(l) +    &    ! coef +1 corresponds to -lwout included in calculation of radnet_l
    4016                      (3.0_wp+1.0_wp) * emiss_surf(l) * sigma_sb * t_surf(l) ** 4 +      & 
    4017                      f_shf  * pt1 +                                                     &
    4018                      lambda_surface * t_wall(nzb_wall,l)
    4019 
    4020 !--         denominator of the prognostic equation
    4021             coef_2 = 4.0_wp * emiss_surf(l) * sigma_sb * t_surf(l) ** 3                 &
    4022                          + lambda_surface + f_shf / exn(k)
    4023 
    4024 !--         implicit solution when the surface layer has no heat capacity,
    4025 !--         otherwise use RK3 scheme.
    4026             t_surf_p(l) = ( coef_1 * dt_3d * tsc(2) + c_surface(l) * t_surf(l) ) /      &
    4027                               ( c_surface(l) + coef_2 * dt_3d * tsc(2) )
    4028 
    4029 !--         add RK3 term
    4030             t_surf_p(l) = t_surf_p(l) + dt_3d * tsc(3) * tt_surface_m(l)
     4906             Ueff   = SQRT( u1**2 + v1**2 + w1**2 )
     4907             httc   = surf_usm_v(l)%roughness_wall(m) *                       &
     4908                       ( 11.8_wp + 4.2_wp * Ueff ) - 4.0_wp
     4909             f_shf  = httc
     4910
     4911!--          add LW up so that it can be removed in prognostic equation
     4912             surf_usm_v(l)%rad_net_l(m) = surf_usm_v(l)%rad_in_sw(m)  -        &
     4913                                          surf_usm_v(l)%rad_out_sw(m) +        &
     4914                                          surf_usm_v(l)%rad_in_lw(m)  -        &
     4915                                          surf_usm_v(l)%rad_out_lw(m)
     4916
     4917!--           numerator of the prognostic equation
     4918              coef_1 = surf_usm_v(l)%rad_net_l(m) +                            & ! coef +1 corresponds to -lwout included in calculation of radnet_l
     4919               ( 3.0_wp + 1.0_wp ) * surf_usm_v(l)%emiss_surf(m) * sigma_sb *  &
     4920                                     t_surf_v(l)%t(m) ** 4 +                   & 
     4921                                     f_shf * pt1 +                             &
     4922                                     lambda_surface * t_wall_v(l)%t(nzb_wall,m)
     4923
     4924!--           denominator of the prognostic equation
     4925              coef_2 = 4.0_wp * surf_usm_v(l)%emiss_surf(m) * sigma_sb *       &
     4926                                t_surf_v(l)%t(m) ** 3                          &
     4927                              + lambda_surface + f_shf / exn(k)
     4928
     4929!--           implicit solution when the surface layer has no heat capacity,
     4930!--           otherwise use RK3 scheme.
     4931              t_surf_v_p(l)%t(m) = ( coef_1 * dt_3d * tsc(2) +                 &
     4932                             surf_usm_v(l)%c_surface(m) * t_surf_v(l)%t(m) ) / &
     4933                           ( surf_usm_v(l)%c_surface(m) + coef_2 * dt_3d * tsc(2) )
     4934
     4935!--           add RK3 term
     4936              t_surf_v_p(l)%t(m) = t_surf_v_p(l)%t(m) + dt_3d * tsc(3) *       &
     4937                                surf_usm_v(l)%tt_surface_m(m)
    40314938           
    4032 !--         calculate true tendency
    4033             stend = (t_surf_p(l) - t_surf(l) - dt_3d * tsc(3) * tt_surface_m(l)) / (dt_3d  * tsc(2))
    4034 
    4035 !--         calculate t_surf tendencies for the next Runge-Kutta step
    4036             IF ( timestep_scheme(1:5) == 'runge' )  THEN
    4037                 IF ( intermediate_timestep_count == 1 )  THEN
    4038                     tt_surface_m(l) = stend
    4039                 ELSEIF ( intermediate_timestep_count <                                  &
    4040                          intermediate_timestep_count_max )  THEN
    4041                     tt_surface_m(l) = -9.5625_wp * stend + 5.3125_wp                    &
    4042                                        * tt_surface_m(l)
    4043                 ENDIF
    4044             ENDIF
    4045 
    4046 !--         in case of fast changes in the skin temperature, it is required to
    4047 !--         update the radiative fluxes in order to keep the solution stable
    4048             IF ( ABS( t_surf_p(l) - t_surf(l) ) > 1.0_wp )  THEN
    4049                force_radiation_call_l = .TRUE.
    4050             ENDIF
    4051            
    4052 !--         for horizontal surfaces is pt(nzb_s_inner(j,i),j,i) = pt_surf.
    4053 !--         there is no equivalent surface gridpoint for vertical surfaces.
    4054 !--         pt(k,j,i) is calculated for all directions in diffusion_s
    4055 !--         using surface and wall heat fluxes
    4056             IF ( d == iroof )  THEN
    4057                pt(nzb_s_inner(j,i),j,i) = t_surf_p(l) / exn(k)
    4058             ENDIF
    4059 
    4060 !--         calculate fluxes
    4061 !--         rad_net_l is never used!           
    4062             rad_net_l(l)     = rad_net_l(l) + 3.0_wp * sigma_sb                         &
    4063                                 * t_surf(l)**4 - 4.0_wp * sigma_sb                      &
    4064                                 * t_surf(l)**3 * t_surf_p(l)
    4065             wghf_eb(l)       = lambda_surface * (t_surf_p(l) - t_wall(nzb_wall,l))
    4066 
    4067 !--         ground/wall/roof surface heat flux
    4068             wshf_eb(l)  = - f_shf  * ( pt1 - t_surf_p(l) )
    4069            
    4070 !--         store kinematic surface heat fluxes for utilization in other processes
    4071 !--         diffusion_s, surface_layer_fluxes,...
    4072             IF ( d == iroof )  THEN
    4073 !--             shf is used in diffusion_s and also
    4074 !--             for calculation of surface layer fluxes
    4075 !--             update for horizontal surfaces
    4076                 shf(j,i) = wshf_eb(l) / rho_cp
    4077             ELSE
    4078 !--             surface heat flux for vertical surfaces
    4079 !--             used in diffusion_s
    4080                 wshf(l) = wshf_eb(l) / rho_cp
    4081             ENDIF
     4939!--           calculate true tendency
     4940              stend = ( t_surf_v_p(l)%t(m) - t_surf_v(l)%t(m) - dt_3d * tsc(3) *&
     4941                        surf_usm_v(l)%tt_surface_m(m) ) / ( dt_3d  * tsc(2) )
     4942
     4943!--           calculate t_surf tendencies for the next Runge-Kutta step
     4944              IF ( timestep_scheme(1:5) == 'runge' )  THEN
     4945                 IF ( intermediate_timestep_count == 1 )  THEN
     4946                    surf_usm_v(l)%tt_surface_m(m) = stend
     4947                 ELSEIF ( intermediate_timestep_count <                        &
     4948                          intermediate_timestep_count_max )  THEN
     4949                    surf_usm_v(l)%tt_surface_m(m) = -9.5625_wp * stend +       &
     4950                                     5.3125_wp * surf_usm_h%tt_surface_m(m)
     4951                 ENDIF
     4952              ENDIF
     4953
     4954!--           in case of fast changes in the skin temperature, it is required to
     4955!--           update the radiative fluxes in order to keep the solution stable
     4956              IF ( ABS( t_surf_v_p(l)%t(m) - t_surf_v(l)%t(m) ) > 1.0_wp )  THEN
     4957                 force_radiation_call_l = .TRUE.
     4958              ENDIF
     4959
     4960!--           calculate fluxes
     4961!--           rad_net_l is never used!           
     4962              surf_usm_v(l)%rad_net_l(m) = surf_usm_v(l)%rad_net_l(m) +        &
     4963                                     3.0_wp * sigma_sb *                       &
     4964                                     t_surf_v(l)%t(m)**4 - 4.0_wp * sigma_sb * &
     4965                                     t_surf_v(l)%t(m)**3 * t_surf_v_p(l)%t(m)
     4966
     4967              surf_usm_v(l)%wghf_eb(m)   = lambda_surface *                    &
     4968                                     ( t_surf_v_p(l)%t(m) - t_wall_v(l)%t(nzb_wall,m) )
     4969
     4970!--           ground/wall/roof surface heat flux
     4971              surf_usm_v(l)%wshf_eb(m)   = - f_shf  * ( pt1 - t_surf_v_p(l)%t(m) )
     4972
     4973!           
     4974!--           store kinematic surface heat fluxes for utilization in other processes
     4975!--           diffusion_s, surface_layer_fluxes,...
     4976              surf_usm_v(l)%shf(m) = surf_usm_v(l)%wshf_eb(m) / rho_cp
     4977
     4978           ENDDO
    40824979
    40834980        ENDDO
    4084        
    4085        
     4981!
     4982!--     Add-up anthropogenic heat, for now only at upward-facing surfaces
    40864983        IF ( usm_anthropogenic_heat  .AND.  &
    40874984             intermediate_timestep_count == intermediate_timestep_count_max )  THEN
    4088 !--         application of the additional anthropogenic heat sources
    4089 !--         we considere the traffic for now so all heat is absorbed
    4090 !--         to the first layer, generalization would be worth
     4985!--        application of the additional anthropogenic heat sources
     4986!--        we considere the traffic for now so all heat is absorbed
     4987!--        to the first layer, generalization would be worth.
    40914988           
    4092 !--         calculation of actual profile coefficient
    4093 !--         ??? check time_since_reference_point ???
    4094             dtime = mod(simulated_time + time_utc_init, 24.0_wp*3600.0_wp)
    4095             dhour = INT(dtime/3600.0_wp)
    4096 !--         linear interpolation of coeficient
    4097             acoef = (REAL(dhour+1,wp)-dtime/3600.0_wp)*aheatprof(dhour) + (dtime/3600.0_wp-REAL(dhour,wp))*aheatprof(dhour+1)
    4098             DO i = nxl, nxr
    4099                 DO j = nys, nyn
    4100                     IF ( aheat(j,i) > 0.0_wp )  THEN
    4101 !--                     TODO the increase of pt in box i,j,nzb_s_inner(j,i)+1 in time dt_3d
    4102 !--                     given to anthropogenic heat aheat*acoef (W*m-2)
    4103 !--                     k = nzb_s_inner(j,i)+1
    4104 !--                     pt(k,j,i) = pt(k,j,i) + aheat(j,i)*acoef*dt_3d/(exn(k)*rho_cp*dz)
    4105 !--                     Instead of this, we can adjust shf in case AH only at surface
    4106                         shf(j,i) = shf(j,i) + aheat(j,i)*acoef * ddx * ddy / rho_cp
    4107                     ENDIF
    4108                 ENDDO
    4109             ENDDO
     4989!--        calculation of actual profile coefficient
     4990!--        ??? check time_since_reference_point ???
     4991           dtime = mod(simulated_time + time_utc_init, 24.0_wp*3600.0_wp)
     4992           dhour = INT(dtime/3600.0_wp)
     4993!--        linear interpolation of coeficient
     4994           acoef = (REAL(dhour+1,wp)-dtime/3600.0_wp)*aheatprof(dhour) + (dtime/3600.0_wp-REAL(dhour,wp))*aheatprof(dhour+1)
     4995
     4996           DO m = 1, surf_usm_h%ns
     4997!
     4998!--           Get indices of respective grid point
     4999              i = surf_usm_h%i(m)
     5000              j = surf_usm_h%j(m)
     5001              k = surf_usm_h%k(m)
     5002
     5003              IF ( aheat(j,i) > 0.0_wp )  THEN
     5004!--              TODO the increase of pt in box i,j,nzb_s_inner(j,i)+1 in time dt_3d
     5005!--              given to anthropogenic heat aheat*acoef (W*m-2)
     5006!--              k = nzb_s_inner(j,i)+1
     5007!--              pt(k,j,i) = pt(k,j,i) + aheat(j,i)*acoef*dt_3d/(exn(k)*rho_cp*dz)
     5008!--              Instead of this, we can adjust shf in case AH only at surface
     5009                 surf_usm_h%shf(m) = surf_usm_h%shf(m) +                       &
     5010                                   aheat(j,i) * acoef * ddx * ddy / rho_cp
     5011              ENDIF
     5012           ENDDO
     5013
    41105014        ENDIF
    41115015       
     
    41135017!--     get the borders from neighbours
    41145018        CALL exchange_horiz( pt, nbgp )
    4115         CALL exchange_horiz_2d( shf )
    4116 
    4117 
    4118 !--    calculation of force_radiation_call:
    4119 !--    Make logical OR for all processes.
    4120 !--    Force radiation call if at least one processor forces it.
    4121        IF ( intermediate_timestep_count == intermediate_timestep_count_max-1 )          &
    4122        THEN
     5019
     5020
     5021!--     calculation of force_radiation_call:
     5022!--     Make logical OR for all processes.
     5023!--     Force radiation call if at least one processor forces it.
     5024        IF ( intermediate_timestep_count == intermediate_timestep_count_max-1 )&
     5025        THEN
    41235026#if defined( __parallel )
    41245027          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    4125               CALL mpi_allreduce( force_radiation_call_l, force_radiation_call,         &
    4126                                   1, MPI_LOGICAL, MPI_LOR, comm2d, ierr )
     5028          CALL MPI_ALLREDUCE( force_radiation_call_l, force_radiation_call,    &
     5029                              1, MPI_LOGICAL, MPI_LOR, comm2d, ierr )
    41275030#else
    41285031          force_radiation_call = force_radiation_call_l
     
    41305033          force_radiation_call_l = .FALSE.
    41315034       ENDIF
    4132 
    41335035
    41345036    END SUBROUTINE usm_surface_energy_balance
     
    41495051     
    41505052#if defined( __nopointer )
    4151        t_surf    = t_surf_p
    4152        t_wall    = t_wall_p
     5053       t_surf_h    = t_surf_h_p
     5054       t_wall_h    = t_wall_h_p
     5055       t_surf_v    = t_surf_v_p
     5056       t_wall_v    = t_wall_v_p
    41535057#else
    41545058       SELECT CASE ( mod_count )
    41555059          CASE ( 0 )
    4156              t_surf  => t_surf_1; t_surf_p  => t_surf_2
    4157              t_wall     => t_wall_1;    t_wall_p     => t_wall_2
     5060!
     5061!--          Horizontal surfaces
     5062             t_surf_h  => t_surf_h_1; t_surf_h_p  => t_surf_h_2
     5063             t_wall_h     => t_wall_h_1;    t_wall_h_p     => t_wall_h_2
     5064!
     5065!--          Vertical surfaces
     5066             t_surf_v  => t_surf_v_1; t_surf_v_p  => t_surf_v_2
     5067             t_wall_v     => t_wall_v_1;    t_wall_v_p     => t_wall_v_2
    41585068          CASE ( 1 )
    4159              t_surf  => t_surf_2; t_surf_p  => t_surf_1
    4160              t_wall     => t_wall_2;    t_wall_p     => t_wall_1
     5069!
     5070!--          Horizontal surfaces
     5071             t_surf_h  => t_surf_h_2; t_surf_h_p  => t_surf_h_1
     5072             t_wall_h     => t_wall_h_2;    t_wall_h_p     => t_wall_h_1
     5073!
     5074!--          Vertical surfaces
     5075             t_surf_v  => t_surf_v_2; t_surf_v_p  => t_surf_v_1
     5076             t_wall_v     => t_wall_v_2;    t_wall_v_p     => t_wall_v_1
    41615077       END SELECT
    41625078#endif
    41635079       
    41645080    END SUBROUTINE usm_swap_timelevel
    4165 
    4166 
    4167 !------------------------------------------------------------------------------!
    4168 ! Description:
    4169 ! ------------
    4170 !
    4171 !> This function applies the kinematic wall heat fluxes
    4172 !> for walls in four directions for all gridboxes in urban layer.
    4173 !> It is called out from subroutine prognostic_equations.
    4174 !> TODO Compare performance with cycle runnig l=startwall,endwall...
    4175 !------------------------------------------------------------------------------!
    4176     SUBROUTINE usm_wall_heat_flux
    4177    
    4178         IMPLICIT NONE
    4179 
    4180         INTEGER(iwp)              ::  i,j,k,d,l             !< running indices
    4181        
    4182         DO l = startenergy, endenergy
    4183             j = surfl(iy,l)
    4184             i = surfl(ix,l)
    4185             k = surfl(iz,l)
    4186             d = surfl(id,l)
    4187             tend(k,j,i) = tend(k,j,i) + wshf(l) * ddxy2(d)
    4188         ENDDO
    4189 
    4190     END SUBROUTINE usm_wall_heat_flux
    4191  
    4192  
    4193 !------------------------------------------------------------------------------!
    4194 ! Description:
    4195 ! ------------
    4196 !
    4197 !> This function applies the kinematic wall heat fluxes
    4198 !> for walls in four directions around the gridbox i,j.
    4199 !> It is called out from subroutine prognostic_equations.
    4200 !------------------------------------------------------------------------------!
    4201     SUBROUTINE usm_wall_heat_flux_ij(i,j)
    4202    
    4203         IMPLICIT NONE
    4204 
    4205         INTEGER(iwp), INTENT(in)  ::  i,j                   !< indices of grid box
    4206         INTEGER(iwp)              ::  ii,jj,k,d,l
    4207        
    4208         DO l = startenergy, endenergy
    4209             jj = surfl(iy,l)
    4210             ii = surfl(ix,l)
    4211             IF ( ii == i  .AND.  jj == j ) THEN
    4212                k = surfl(iz,l)
    4213                IF ( k >=  nzb_s_inner(j,i)+1  .AND.  k <=  nzb_s_outer(j,i) ) THEN
    4214                   d = surfl(id,l)
    4215                   IF ( d >= 1 .and. d <= 4 )   THEN
    4216                      tend(k,j,i) = tend(k,j,i) + wshf(l) * ddxy2(d)
    4217                   ENDIF
    4218                ENDIF
    4219             ENDIF
    4220         ENDDO
    4221 
    4222     END SUBROUTINE usm_wall_heat_flux_ij
    4223  
    42245081
    42255082!------------------------------------------------------------------------------!
     
    42425099       DO  i = 0, io_blocks-1
    42435100          IF ( i == io_group )  THEN
    4244              WRITE ( 14 )  't_surf                        '
     5101
     5102             WRITE ( 14 )  't_surf_h                      '
    42455103#if defined( __nopointer )             
    4246              WRITE ( 14 )  t_surf
     5104             WRITE ( 14 )  t_surf_h
    42475105#else
    4248              WRITE ( 14 )  t_surf_1
     5106             WRITE ( 14 )  t_surf_h_1
    42495107#endif
    4250              WRITE ( 14 )  't_wall                        '
     5108             WRITE ( 14 )  't_surf_v(0)                   '
    42515109#if defined( __nopointer )             
    4252              WRITE ( 14 )  t_wall
     5110             WRITE ( 14 )  t_surf_v(0)%t
    42535111#else
    4254              WRITE ( 14 )  t_wall_1
     5112             WRITE ( 14 )  t_surf_v_1(0)%t
     5113#endif
     5114             WRITE ( 14 )  't_surf_v(1)                   '
     5115#if defined( __nopointer )             
     5116             WRITE ( 14 )  t_surf_v(1)%t
     5117#else
     5118             WRITE ( 14 )  t_surf_v_1(1)%t
     5119#endif
     5120             WRITE ( 14 )  't_surf_v(2)                   '
     5121#if defined( __nopointer )             
     5122             WRITE ( 14 )  t_surf_v(2)%t
     5123#else
     5124             WRITE ( 14 )  t_surf_v_1(2)%t
     5125#endif
     5126             WRITE ( 14 )  't_surf_v(3)                   '
     5127#if defined( __nopointer )             
     5128             WRITE ( 14 )  t_surf_v(3)%t
     5129#else
     5130             WRITE ( 14 )  t_surf_v_1(3)%t
     5131#endif
     5132             WRITE ( 14 )  't_wall_h                      '
     5133#if defined( __nopointer )             
     5134             WRITE ( 14 )  t_wall_h
     5135#else
     5136             WRITE ( 14 )  t_wall_h_1
     5137#endif
     5138             WRITE ( 14 )  't_wall_v(0)                   '
     5139#if defined( __nopointer )             
     5140             WRITE ( 14 )  t_wall_v(0)%t
     5141#else
     5142             WRITE ( 14 )  t_wall_v_1(0)%t
     5143#endif
     5144             WRITE ( 14 )  't_wall_v(1)                   '
     5145#if defined( __nopointer )             
     5146             WRITE ( 14 )  t_wall_v(1)%t
     5147#else
     5148             WRITE ( 14 )  t_wall_v_1(1)%t
     5149#endif
     5150             WRITE ( 14 )  't_wall_v(2)                   '
     5151#if defined( __nopointer )             
     5152             WRITE ( 14 )  t_wall_v(2)%t
     5153#else
     5154             WRITE ( 14 )  t_wall_v_1(2)%t
     5155#endif
     5156             WRITE ( 14 )  't_wall_v(3)                   '
     5157#if defined( __nopointer )             
     5158             WRITE ( 14 )  t_wall_v(3)%t
     5159#else
     5160             WRITE ( 14 )  t_wall_v_1(3)%t
    42555161#endif
    42565162             WRITE ( 14 )  '*** end usm ***               '
Note: See TracChangeset for help on using the changeset viewer.