Changeset 4159 for palm/trunk/SOURCE/init_grid.f90
- Timestamp:
- Aug 15, 2019 1:31:35 PM (5 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/init_grid.f90 (modified) (12 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/init_grid.f90
r4144 r4159 25 25 ! ----------------- 26 26 ! $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 27 31 ! relational operators .EQ., .NE., etc. replaced by ==, /=, etc. 28 32 ! … … 1408 1412 !-- need to be remove to avoid steep edges at the child-domain boundaries. 1409 1413 IF ( input_pids_static ) THEN 1414 1410 1415 #if defined( __parallel ) 1411 1416 CALL MPI_ALLREDUCE( MINVAL( terrain_height_f%var ), oro_min, 1, & … … 1414 1419 oro_min = MINVAL( terrain_height_f%var ) 1415 1420 #endif 1416 1417 1421 terrain_height_f%var = terrain_height_f%var - oro_min 1418 1422 ! … … 1575 1579 1576 1580 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) == & 1580 1585 build_ids_final(nr) ) ) 1581 1586 ENDDO … … 1583 1588 #if defined( __parallel ) 1584 1589 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,& 1586 1591 MPI_MAX, comm2d, ierr ) 1587 1592 ENDIF … … 1591 1596 ! 1592 1597 !-- 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 1594 1602 oro_max_l = 0.0 1595 1603 DO nr = 1, SIZE(build_ids_final) … … 1625 1633 !-- attributes will not be correct as given surface information 1626 1634 !-- 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. 1629 1636 IF ( zu(k) - ocean_offset <= terrain_height_f%var(j,i) ) THEN 1630 1637 topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 ) … … 1636 1643 !-- 3D buildings require separate treatment. 1637 1644 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 1641 1656 topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 ) 1642 1657 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 )1646 1658 ENDIF 1647 1659 ENDIF 1648 1660 ENDIF 1649 1661 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 1650 1682 ! 1651 1683 !-- Map 3D buildings onto terrain height. … … 1669 1701 IF ( building_type_f%var(j,i) /= 7 ) THEN 1670 1702 DO k = topo_top_index + 1, nzt + 1 1671 IF ( z w(k) - ocean_offset <= oro_max(nr) ) THEN1703 IF ( zu(k) - ocean_offset <= oro_max(nr) ) THEN 1672 1704 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 ) 1674 1706 ENDIF 1675 1707 ENDDO … … 1678 1710 !-- lower start index where building starts. 1679 1711 DO k = nzb, nzt 1680 IF ( z w(k) - ocean_offset <= oro_max(nr) ) &1712 IF ( zu(k) - ocean_offset <= oro_max(nr) ) & 1681 1713 topo_top_index = k 1682 1714 ENDDO … … 1690 1722 IF ( buildings_f%var_3d(k2,j,i) == 1 ) THEN 1691 1723 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 )1693 1724 topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 2 ) 1694 1725 ENDIF … … 1709 1740 ELSE 1710 1741 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. 1711 1749 topo_3d = IBSET( topo_3d, 0 ) 1712 1750 topo_3d(nzb,:,:) = IBCLR( topo_3d(nzb,:,:), 0 ) 1751 topo_3d(nzb,:,:) = IBSET( topo_3d(nzb,:,:), 1 ) 1713 1752 DO i = nxl, nxr 1714 1753 DO j = nys, nyn
Note: See TracChangeset
for help on using the changeset viewer.
