Ignore:
Timestamp:
Aug 15, 2019 1:31:35 PM (22 months ago)
Author:
suehring
Message:

Revision of topography processing to have a consistent treatment of 2D and 3D buildings (init_grid, surface_mod); Bugfix in indoor model in case of non grid-resolved buildings

File:
1 edited

Legend:

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

    r4144 r4159  
    2525! -----------------
    2626! $Id$
     27! Revision of topography processing. This was not consistent between 2D and 3D
     28! buildings.
     29!
     30! 4144 2019-08-06 09:11:47Z raasch
    2731! relational operators .EQ., .NE., etc. replaced by ==, /=, etc.
    2832!
     
    14081412!-- need to be remove to avoid steep edges at the child-domain boundaries.
    14091413    IF ( input_pids_static )  THEN
     1414   
    14101415#if defined( __parallel )
    14111416       CALL MPI_ALLREDUCE( MINVAL( terrain_height_f%var ), oro_min, 1,         &
     
    14141419       oro_min = MINVAL( terrain_height_f%var )
    14151420#endif
    1416 
    14171421       terrain_height_f%var = terrain_height_f%var - oro_min
    14181422!                           
     
    15751579
    15761580          DO  nr = 1, SIZE(build_ids_final)
    1577              oro_max_l(nr) = MAXVAL(                                              &
    1578                               MERGE( terrain_height_f%var, 0.0_wp,                &
    1579                                      building_id_f%var(nys:nyn,nxl:nxr) ==      &
     1581             oro_max_l(nr) = MAXVAL(                                           &
     1582                              MERGE( terrain_height_f%var(nys:nyn,nxl:nxr),    &
     1583                                     0.0_wp,                                   &
     1584                                     building_id_f%var(nys:nyn,nxl:nxr) ==     &
    15801585                                     build_ids_final(nr) ) )
    15811586          ENDDO
     
    15831588#if defined( __parallel )   
    15841589          IF ( SIZE(build_ids_final) >= 1 ) THEN
    1585              CALL MPI_ALLREDUCE( oro_max_l, oro_max, SIZE( oro_max ), MPI_REAL,   &
     1590             CALL MPI_ALLREDUCE( oro_max_l, oro_max, SIZE( oro_max ), MPI_REAL,&
    15861591                                 MPI_MAX, comm2d, ierr )
    15871592          ENDIF
     
    15911596!
    15921597!--       Finally, determine discrete grid height of maximum orography occupied
    1593 !--       by a building. Use all-or-nothing approach, i.e. a grid box is either
     1598!--       by a building. Use all-or-nothing approach, i.e. if terrain
     1599!--       exceeds the scalar level the grid box is fully terrain and the
     1600!--       maximum terrain is set to the zw level.
     1601!--       terrain or
    15941602          oro_max_l = 0.0
    15951603          DO  nr = 1, SIZE(build_ids_final)
     
    16251633!--             attributes will not be correct as given surface information
    16261634!--             will not be in accordance to the classified grid points.
    1627 !--             Hence, in this case, de-flag the grid point and give it
    1628 !--             urban type instead.
     1635!--             Hence, in this case, also a building flag.
    16291636                IF ( zu(k) - ocean_offset <= terrain_height_f%var(j,i) )  THEN
    16301637                    topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
     
    16361643!--             3D buildings require separate treatment.
    16371644                IF ( buildings_f%from_file  .AND.  buildings_f%lod == 1 )  THEN
    1638                    IF ( building_id_f%var(j,i) /= building_id_f%fill )  THEN       
    1639                       IF ( zu(k) - ocean_offset <=                             &
    1640                            oro_max(nr) + buildings_f%var_2d(j,i) )  THEN
     1645!
     1646!--                Fill-up the terrain to the level of maximum orography
     1647!--                within the building-covered area.
     1648                   IF ( building_id_f%var(j,i) /= building_id_f%fill )  THEN
     1649!
     1650!--                   Note, oro_max is always on zw level                   
     1651                      IF ( zu(k) - ocean_offset < oro_max(nr) )  THEN
     1652                         topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
     1653                         topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 1 )
     1654                      ELSEIF ( zu(k) - ocean_offset <=                         &
     1655                               oro_max(nr) + buildings_f%var_2d(j,i) )  THEN
    16411656                         topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
    16421657                         topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 2 )
    1643 !
    1644 !--                      De-flag grid point of type natural. See comment above.
    1645                          topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 1 )
    16461658                      ENDIF
    16471659                   ENDIF
    16481660                ENDIF
    16491661             ENDDO
     1662!
     1663!--          Special treatment for non grid-resolved buildings. This case,
     1664!--          the uppermost terrain grid point is flagged as building as well
     1665!--          well, even though no building exists at all. However, the
     1666!--          surface element will be identified as urban-surface and the
     1667!--          input data provided by the drivers is consistent to the surface
     1668!--          classification. Else, all non grid-resolved buildings would vanish
     1669!--          and identified as terrain grid points, which, however, won't be
     1670!--          consistent with the input data.
     1671             IF ( buildings_f%from_file  .AND.  buildings_f%lod == 1 )  THEN
     1672                IF ( building_id_f%var(j,i) /= building_id_f%fill )  THEN
     1673                   DO  k = nzb, nzt
     1674                      IF( zw(k) - ocean_offset == oro_max(nr) )  THEN
     1675                         IF ( buildings_f%var_2d(j,i) <= zu(k+1) - zw(k) )  THEN
     1676                            topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 2 )
     1677                         ENDIF
     1678                      ENDIF
     1679                   ENDDO
     1680                ENDIF
     1681             ENDIF
    16501682!
    16511683!--          Map 3D buildings onto terrain height. 
     
    16691701                      IF ( building_type_f%var(j,i) /= 7 )  THEN
    16701702                         DO k = topo_top_index + 1, nzt + 1     
    1671                             IF ( zw(k) - ocean_offset <= oro_max(nr) )  THEN
     1703                            IF ( zu(k) - ocean_offset <= oro_max(nr) )  THEN
    16721704                               topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
    1673                                topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 2 )
     1705                               topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 1 )
    16741706                            ENDIF
    16751707                         ENDDO       
     
    16781710!--                      lower start index where building starts.
    16791711                         DO  k = nzb, nzt
    1680                             IF ( zw(k) - ocean_offset <= oro_max(nr) )         &
     1712                            IF ( zu(k) - ocean_offset <= oro_max(nr) )         &
    16811713                               topo_top_index = k
    16821714                         ENDDO
     
    16901722                         IF ( buildings_f%var_3d(k2,j,i) == 1 )  THEN
    16911723                            topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
    1692                             topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 1 )
    16931724                            topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 2 )
    16941725                         ENDIF
     
    17091740    ELSE
    17101741       ocean_offset     = MERGE( zw(0), 0.0_wp, ocean_mode )
     1742!
     1743!--    Initialize topography bit 0 (indicates obstacle) everywhere to zero
     1744!--    and clear all grid points at nzb, where alway a surface is defined.
     1745!--    Further, set also bit 1 (indicates terrain) at nzb, which is further
     1746!--    used for masked data output and further processing. Note, in the
     1747!--    ASCII case no distinction is made between buildings and terrain,
     1748!--    so that setting of bit 1 and 2 at the same time has no effect.
    17111749       topo_3d          = IBSET( topo_3d, 0 )
    17121750       topo_3d(nzb,:,:) = IBCLR( topo_3d(nzb,:,:), 0 )
     1751       topo_3d(nzb,:,:) = IBSET( topo_3d(nzb,:,:), 1 )
    17131752       DO  i = nxl, nxr
    17141753          DO  j = nys, nyn
Note: See TracChangeset for help on using the changeset viewer.