Topography implementation

Topography masking

With revision -r2232, the topography implementation is completely revised. Starting from this revision, overhanging structures (e.g. bridges, ceilings or tunnels) are allowed, i.e. topography does not necessarily be surface-mounted.

The topography implementation described in Sect. boundary conditions allows the use of 2-D topography height data (if no overhanging structures should be considered), or 3-D topography information. The topography input data has to be provided within a rastered ASCII or NetCDF file. Rastered 3-D topography consists of values either 1 (within topography) or 0 (outside topography), and has to be provided up to the level of the highest present topography structures. After reading and mapping of topography data to the 3-D grid in PALM, flag arrays are set in order to mask grid boxes within the integration loops accordingly.

The prognostic terms are executed overall, i.e. also within topography. Masking of topography is done via Fortran integer bit flags, where a certain bit position indicates whether the respective grid point belongs to an obstacle (bit is 0) or to the atmosphere (bit is 1). Different bit positions indicate different masking on the staggered grid, e.g. on the u, v, w grid, equivalent to the former nzb_u/v/w_inner arrays. Further, regions where special surface-bounded code is executed, is also masked via special bit flags. A list of the respective bit flags and their meaning concerning the deprecated implementation via the nzb_u/v/w_inner arrays can be found in the following Table. Bit flags are stored in array wall_flags_0.

Bit position	Meaning
0	mask topography on scalar-grid ( former nzb_s_inner )
1	mask topography on u-grid ( former nzb_u_inner )
2	mask topography on v-grid ( former nzb_v_inner )
3	mask topography on w-grid ( former nzb_w_inner )
8	mask regions where surface-bounded code for u, v, w, and scalars is executed ( combines information of former nzb_s_outer, nzb_u_outer, nzb_v_outer and nzb_w_outer arrays )
9	mask regions where model-top fluxes are applied ( former nzt_diff )
12	flags upward-facing surfaces on scalar-grid
13	flags downward-facing surfaces on scalar-grid
14	flags upward-facing surfaces on u-grid
15	flags downward-facing surfaces on u-grid
17	flags upward-facing surfaces on v-grid
17	flags downward-facing surfaces on v-grid
18	flags upward-facing surfaces on w-grid
19	flags downward-facing surfaces on w-grid
20	mask topography on u-grid and one grid point above, used only for adding random perturbations ( former nzb_u_inner+1 )
21	mask topography on v-grid and one grid point above, used only for adding random perturbations ( former nzb_v_inner+1 )
24	mask regions where surface-bounded code is executed on scalar-grid ( former nzb_s_outer )
25	mask regions where surface-bounded code is executed on scalar-grid ( former nzb_s_outer+1 )
26	mask regions where surface-bounded code is executed on u-grid ( former nzb_u_outer )
27	mask regions where surface-bounded code is executed on v-grid ( former nzb_v_outer )
28	mask regions where surface-bounded code is executed on w-grid ( former nzb_w_outer )
29	mask regions where surface-bounded code is executed on s-grid ( former nzb_diff_s_outer - 1 )
30	mask regions where surface-bounded code is executed on s-grid ( former nzb_diff_s_outer )

The following example illustrates the general realization of topography masking in the code using Bit flags:

Here, the intrinsic Fortran function MERGE returns 1 if wall_flags_0(k,j,i) at bit position zero is 1 (BTEST returns .TRUE.), else, MERGE would return 0 and no term is effectively added to the tend array. A similar approach is also realized for the 5th-order advection-scheme to degrade the order near surfaces.

With revision -r 2232, it is also possible to prescribe topography via NetCDF data format (ASCII still exists). Topography input via NetCDF enables two different modes.
1) Topography height can be prescribed for each grid point, similar to the rastered ASCII iput.
2) Topography is prescribed via a three-dimension Integer input array which is 1 for topography grid points and 0 everywhere else.
This enables the user to prescribe three-dimensional topography structures like bridges, etc. . The input mode is determined internally according to the NetCDF file attribute lod, which is part of the NetCDF input file. lod = 1 for method 1, lod = 2 for method 2.

Remark concerning still existing 2-D index arrays: Even though 2-D index arrays, as e.g. nzb_s_inner, are removed from the integration loops, they can still be used by the user for a while, e.g. for setting user-defined topography.

Treatment of holes

In case there are holes in the topography that are resolved by only one grid point, the topography array in filtered so that such holes are filled up to the minimum topography height of the directly adjoining grid points in the x- and the y-direction. This is necessary, because for such chimney-like features resolved by one grid point, the continuity equation is not fulfilled on a discrete grid, as the only degree of freedom for the pressure solver is the vertical. Such holes are suspected to lead to unrealistic velocity blow-ups, hence, they are filled up.

Surface-bounded code

To execute surface-bounded code, a new Fortran data-structure surf_type is introduced, encompassing all required surface attributes. Arrays belonging to these data-structure are allocated for the exact number of surface elements on a local processor, and loops can run on the exact number of surface elements without any need of IF-ELSE clauses. Following example shows the general structure of the data-structure:

Using this data structure approach, new arrays for different surface types are defined:

The term def, lsm, usm indicates the type of surface, i.e. default-type, natural-type, or urban-type, respectively, while the term at the end indicates either horizontal h or vertical surfaces v. For the different kinds of surfaces, partly different code is executed. For example, for natural- and urban-type surfaces, an energy-balance solver is applied according to the land-surface and urban-surface model, respectively, in order to model surface fluxes of sensible and latent heat.

Further, in addition to the different types of surfaces, surface elements are further distinguished according to its orientation:
surf_X_h(0) encompass all horizontally upward-facing surfaces of its respective type,
surf_X_h(1) encompass all horizontally downward-facing surfaces of its respective type,
surf_X_h(2) encompass the model top grid points.
Please note, downward-facing (and model top) surfaces belong always to the default type at the moment.

Furthermore,
surf_X_v(0) encompass all vertically northward-facing surfaces of its respective type,
surf_X_v(1) encompass all vertically southward-facing surfaces of its respective type,
surf_X_v(2) encompass all vertically eastward-facing surfaces of its respective type,
surf_X_v(3) encompass all vertically westward-facing surfaces of its respective type.

Remark: At the moment it is only possible to have one surface type at once, .e.g. only natural surfaces treated by the land-surface model. In the future, however, different surface types can co-exist next to each other, i.e. it is possible to define areas where the urban-surface model and areas where the land-surface model is executed at the same time.

Outline of technical realization

At first, the number of surface elements of its respective type and orientation is counted and stored within the data-structure:

In the following, surface attributes are allocated for the exact number of surface elements of the respective type:

Here, i, j, k are the corresponding indices linking the surface element to the 3D grid, and shf is the surface heat flux. Please note, there are many other attributes allocated.

A certain variable within the data-structure can be accessed:

Finally, the resulting fluxes are added to the prognostic terms in following way:

Summarized, the surface data-structure surf_type is used to store and access all surface related quantities in an optimized memory demand. Several adjoining surfaces at a grid cell (e.g. at corners) are treated separately. The data-structure is also applied in the flat case without any topography.

Remarks

In order to modify surface properties at given grid point (j,i), it is recommended to follow the subsequent code structure to obtain the type of the respective surface at given (j,i) and the corresponding index of the surface element in the data type.

Here, it it assumed that no overhanging structures exit. For vertical surfaces or in case of overhanging structures, the index m must be looped from start_index(j,i) to end_index(j,i).