source: palm/trunk/SOURCE/init_grid.f90 @ 4420

Last change on this file since 4420 was 4414, checked in by suehring, 4 years ago

Remove deprecated topography arrays; Move basic initialization of numerics into an extra module interface

  • Property svn:keywords set to Id
File size: 117.6 KB
Line 
1!> @file init_grid.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2020 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: init_grid.f90 4414 2020-02-19 20:16:04Z maronga $
27! - Remove deprecated topography arrays nzb_s_inner, nzb_u_inner, etc.
28! - Move initialization of boundary conditions and multigrid into an extra
29!   module interface.
30!
31! 4386 2020-01-27 15:07:30Z Giersch
32! Allocation statements, comments, naming of variables revised and _wp added to
33! real type values
34!
35! 4360 2020-01-07 11:25:50Z suehring
36! Revise error messages for generic tunnel setup.
37!
38! 4346 2019-12-18 11:55:56Z motisi
39! Introduction of wall_flags_total_0, which currently sets bits based on static
40! topography information used in wall_flags_static_0
41!
42! 4340 2019-12-16 08:17:03Z Giersch
43! Topography closed channel flow with symmetric boundaries implemented
44!
45! 4329 2019-12-10 15:46:36Z motisi
46! Renamed wall_flags_0 to wall_flags_static_0
47!
48! 4328 2019-12-09 18:53:04Z suehring
49! Minor change in nzb_max computation. Commentation added.
50!
51! 4314 2019-11-29 10:29:20Z suehring
52! Set additional topography flag 4 to mark topography grid points emerged
53! from the filtering process.
54!
55! 4294 2019-11-13 18:34:16Z suehring
56! Bugfix, always set bit 5 and 6 of wall_flags, indicating terrain- and
57! building surfaces in all  cases, in order to enable terrain-following output
58! also when no land- or urban-surface model is applied.
59!
60! 4265 2019-10-15 16:16:24Z suehring
61! Bugfix for last commit, exchange oro_max variable only when it is allocated
62! (not necessarily the case when topography is input from ASCII file).
63!
64! 4245 2019-09-30 08:40:37Z pavelkrc
65! Store oro_max (building z-offset) in 2D for building surfaces
66!
67! 4189 2019-08-26 16:19:38Z suehring
68! - Add check for proper setting of namelist parameter topography
69! - Set flag to indicate land surfaces in case no topography is provided
70!
71! 4182 2019-08-22 15:20:23Z scharf
72! Corrected "Former revisions" section
73!
74! 4168 2019-08-16 13:50:17Z suehring
75! Pre-calculate topography top index and store it on an array (replaces former
76! functions get_topography_top_index)
77!
78! 4159 2019-08-15 13:31:35Z suehring
79! Revision of topography processing. This was not consistent between 2D and 3D
80! buildings.
81!
82! 4144 2019-08-06 09:11:47Z raasch
83! relational operators .EQ., .NE., etc. replaced by ==, /=, etc.
84!
85! 4115 2019-07-24 12:50:49Z suehring
86! Bugfix in setting near-surface flag 24, inidicating wall-bounded grid points
87!
88! 4110 2019-07-22 17:05:21Z suehring
89! - Separate initialization of advection flags for momentum and scalars.
90! - Change subroutine interface for ws_init_flags_scalar to pass boundary flags
91!
92! 4109 2019-07-22 17:00:34Z suehring
93! Fix bad commit
94!
95! 3926 2019-04-23 12:56:42Z suehring
96! Minor bugfix in building mapping when all building IDs in the model domain
97! are missing
98!
99! 3857 2019-04-03 13:00:16Z knoop
100! In projection of non-building 3D objects onto numerical grid remove
101! dependency on building_type
102!
103! 3763 2019-02-25 17:33:49Z suehring
104! Replace work-around for ghost point exchange of 1-byte arrays with specific
105! routine as already done in other routines
106!
107! 3761 2019-02-25 15:31:42Z raasch
108! unused variables removed
109!
110! 3661 2019-01-08 18:22:50Z suehring
111! Remove setting of nzb_max to nzt at non-cyclic boundary PEs, instead,
112! order degradation of advection scheme is handeled directly in advec_ws
113!
114! 3655 2019-01-07 16:51:22Z knoop
115! Comment added
116!
117! Revision 1.1  1997/08/11 06:17:45  raasch
118! Initial revision (Testversion)
119!
120!
121! Description:
122! -----------------------------------------------------------------------------!
123!> Creating grid depending constants
124!> @todo: Rearrange topo flag list
125!> @todo: reference 3D buildings on top of orography is not tested and may need
126!>        further improvement for steep slopes
127!> @todo: Use more advanced setting of building type at filled holes
128!------------------------------------------------------------------------------!
129 SUBROUTINE init_grid
130
131    USE arrays_3d,                                                             &
132        ONLY:  dd2zu, ddzu, ddzu_pres, ddzw, dzu, dzw, x, xu, y, yv, zu, zw
133
134    USE control_parameters,                                                    &
135        ONLY:  constant_flux_layer, dz, dz_max, dz_stretch_factor,             &
136               dz_stretch_factor_array, dz_stretch_level, dz_stretch_level_end,&
137               dz_stretch_level_end_index, dz_stretch_level_start_index,       &
138               dz_stretch_level_start, ibc_uv_b, message_string,               &
139               number_stretch_level_end,                                       &
140               number_stretch_level_start,                                     &
141               ocean_mode,                                                     &
142               psolver,                                                        &
143               symmetry_flag,                                                  &
144               topography,                                                     &
145               use_surface_fluxes
146
147    USE grid_variables,                                                        &
148        ONLY:  ddx, ddx2, ddy, ddy2, dx, dx2, dy, dy2, zu_s_inner, zw_w_inner
149
150    USE indices,                                                               &
151        ONLY:  nbgp,                                                           &
152               nx,                                                             &
153               nxl,                                                            &
154               nxlg,                                                           &
155               nxr,                                                            &
156               nxrg,                                                           &
157               ny,                                                             &
158               nyn,                                                            &
159               nyng,                                                           &
160               nys,                                                            &
161               nysg,                                                           &
162               nz,                                                             &
163               nzb,                                                            &
164               nzb_diff,                                                       &
165               nzb_max,                                                        &
166               nzt,                                                            &
167               topo_top_ind,                                                   &
168               topo_min_level
169
170    USE kinds
171
172    USE pegrid
173
174    USE vertical_nesting_mod,                                                  &
175        ONLY:  vnested, vnest_init_grid
176
177    IMPLICIT NONE
178
179    INTEGER(iwp) ::  i             !< index variable along x
180    INTEGER(iwp) ::  j             !< index variable along y
181    INTEGER(iwp) ::  k             !< index variable along z
182    INTEGER(iwp) ::  k_top         !< topography top index on local PE
183    INTEGER(iwp) ::  n             !< loop variable for stretching
184    INTEGER(iwp) ::  number_dz     !< number of user-specified dz values       
185    INTEGER(iwp) ::  nzb_local_max !< vertical grid index of maximum topography height
186    INTEGER(iwp) ::  nzb_local_min !< vertical grid index of minimum topography height
187
188    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE ::  topo !< input array for 3D topography and dummy array for setting "outer"-flags
189
190    REAL(wp) ::  dz_level_end  !< distance between calculated height level for u/v-grid and user-specified end level for stretching
191    REAL(wp) ::  dz_stretched  !< stretched vertical grid spacing
192   
193    REAL(wp), DIMENSION(:), ALLOCATABLE ::  min_dz_stretch_level_end !< Array that contains all minimum heights where the stretching can end
194
195
196!
197!-- Calculation of horizontal array bounds including ghost layers
198    nxlg = nxl - nbgp
199    nxrg = nxr + nbgp
200    nysg = nys - nbgp
201    nyng = nyn + nbgp
202
203!
204!-- Allocate grid arrays
205    ALLOCATE( x(0:nx) )
206    ALLOCATE( xu(0:nx) )
207   
208    DO i = 0, nx
209       xu(i) = i * dx
210       x(i)  = i * dx + 0.5_wp * dx
211    ENDDO
212
213    ALLOCATE( y(0:ny) )
214    ALLOCATE( yv(0:ny) )
215   
216    DO j = 0, ny
217       yv(j) = j * dy
218       y(j)  = j * dy + 0.5_wp * dy
219    ENDDO
220
221    ALLOCATE( ddzu(1:nzt+1) )
222    ALLOCATE( ddzw(1:nzt+1) )
223    ALLOCATE( dd2zu(1:nzt) )
224    ALLOCATE( dzu(1:nzt+1) )
225    ALLOCATE( dzw(1:nzt+1) )
226    ALLOCATE( zu(nzb:nzt+1) )
227    ALLOCATE( zw(nzb:nzt+1) )
228
229!
230!-- For constructing an appropriate grid, the vertical grid spacing dz has to
231!-- be specified with a non-negative value in the parameter file
232    IF ( dz(1) == -1.0_wp )  THEN
233       message_string = 'missing dz'
234       CALL message( 'init_grid', 'PA0200', 1, 2, 0, 6, 0 ) 
235    ELSEIF ( dz(1) <= 0.0_wp )  THEN
236       WRITE( message_string, * ) 'dz=',dz(1),' <= 0.0'
237       CALL message( 'init_grid', 'PA0201', 1, 2, 0, 6, 0 )
238    ENDIF
239
240!
241!-- Initialize dz_stretch_level_start with the value of dz_stretch_level
242!-- if it was set by the user
243    IF ( dz_stretch_level /= -9999999.9_wp ) THEN
244       dz_stretch_level_start(1) = dz_stretch_level
245    ENDIF
246       
247!
248!-- Determine number of dz values and stretching levels specified by the
249!-- user to allow right controlling of the stretching mechanism and to
250!-- perform error checks. The additional requirement that dz /= dz_max
251!-- for counting number of user-specified dz values is necessary. Otherwise
252!-- restarts would abort if the old stretching mechanism with dz_stretch_level
253!-- is used (Attention: The user is not allowed to specify a dz value equal
254!-- to the default of dz_max = 999.0).
255    number_dz = COUNT( dz /= -1.0_wp .AND. dz /= dz_max)
256    number_stretch_level_start = COUNT( dz_stretch_level_start /=              &
257                                       -9999999.9_wp )
258    number_stretch_level_end = COUNT( dz_stretch_level_end /=                  &
259                                      9999999.9_wp )
260
261!
262!-- The number of specified end levels +1 has to be the same as the number
263!-- of specified dz values
264    IF ( number_dz /= number_stretch_level_end + 1 ) THEN
265       WRITE( message_string, * ) 'The number of values for dz = ',            &
266                                   number_dz, 'has to be the same as& ',       &
267                                   'the number of values for ',                &
268                                   'dz_stretch_level_end + 1 = ',              &
269                                   number_stretch_level_end+1
270          CALL message( 'init_grid', 'PA0156', 1, 2, 0, 6, 0 )
271    ENDIF
272   
273!
274!-- The number of specified start levels has to be the same or one less than
275!-- the number of specified dz values
276    IF ( number_dz /= number_stretch_level_start + 1 .AND.                     &
277         number_dz /= number_stretch_level_start ) THEN
278       WRITE( message_string, * ) 'The number of values for dz = ',            &
279                                   number_dz, 'has to be the same as or one ', &
280                                   'more than& the number of values for ',     &
281                                   'dz_stretch_level_start = ',                &
282                                   number_stretch_level_start
283          CALL message( 'init_grid', 'PA0211', 1, 2, 0, 6, 0 )
284    ENDIF
285   
286!-- The number of specified start levels has to be the same or one more than
287!-- the number of specified end levels
288    IF ( number_stretch_level_start /= number_stretch_level_end + 1 .AND.      &
289         number_stretch_level_start /= number_stretch_level_end ) THEN
290       WRITE( message_string, * ) 'The number of values for ',                 &
291                                  'dz_stretch_level_start = ',                 &
292                                   dz_stretch_level_start, 'has to be the ',   &
293                                   'same or one more than& the number of ',    &
294                                   'values for dz_stretch_level_end = ',       &
295                                   number_stretch_level_end
296          CALL message( 'init_grid', 'PA0216', 1, 2, 0, 6, 0 )
297    ENDIF
298
299!
300!-- Initialize dz for the free atmosphere with the value of dz_max
301    IF ( dz(number_stretch_level_start+1) == -1.0_wp .AND.                     &
302         number_stretch_level_start /= 0 ) THEN
303       dz(number_stretch_level_start+1) = dz_max
304    ENDIF
305       
306!
307!-- Initialize the stretching factor if (infinitely) stretching in the free
308!-- atmosphere is desired (dz_stretch_level_end was not specified for the
309!-- free atmosphere)
310    IF ( number_stretch_level_start == number_stretch_level_end + 1 ) THEN
311       dz_stretch_factor_array(number_stretch_level_start) =                   &
312       dz_stretch_factor
313    ENDIF
314   
315!
316!-- Allocation of arrays for stretching
317    ALLOCATE( min_dz_stretch_level_end(number_stretch_level_start) )
318
319!
320!-- Define the vertical grid levels. Start with atmosphere branch
321    IF ( .NOT. ocean_mode )  THEN
322   
323!
324!--    The stretching region has to be large enough to allow for a smooth
325!--    transition between two different grid spacings. The number 4 is an
326!--    empirical value
327       DO n = 1, number_stretch_level_start
328          min_dz_stretch_level_end(n) = dz_stretch_level_start(n) +            &
329                                        4 * MAX( dz(n),dz(n+1) )
330       ENDDO
331
332       IF ( ANY( min_dz_stretch_level_end(1:number_stretch_level_start) >      &
333                 dz_stretch_level_end(1:number_stretch_level_start) ) ) THEN
334             message_string= 'Each dz_stretch_level_end has to be larger ' //  &
335                             'than its corresponding value for &' //           &
336                             'dz_stretch_level_start + 4*MAX(dz(n),dz(n+1)) '//&
337                             'to allow for smooth grid stretching'
338             CALL message( 'init_grid', 'PA0224', 1, 2, 0, 6, 0 )
339       ENDIF
340       
341!
342!--    Stretching must not be applied within the surface layer
343!--    (first two grid points). For the default case dz_stretch_level_start
344!--    is negative. Therefore the absolut value is checked here.
345       IF ( ANY( ABS( dz_stretch_level_start ) <= dz(1) * 1.5_wp ) ) THEN
346          WRITE( message_string, * ) 'Each dz_stretch_level_start has to be ',&
347                                     'larger than ', dz(1) * 1.5
348             CALL message( 'init_grid', 'PA0226', 1, 2, 0, 6, 0 )
349       ENDIF
350
351!
352!--    The stretching has to start and end on a grid level. Therefore
353!--    user-specified values are mapped to the next lowest level. The 
354!--    calculation of the first level is realized differently just because of
355!--    historical reasons (the advanced/new stretching mechanism was realized 
356!--    in a way that results don't change if the old parameters
357!--    dz_stretch_level, dz_stretch_factor and dz_max are used)
358       IF ( number_stretch_level_start /= 0 ) THEN
359          dz_stretch_level_start(1) = INT( (dz_stretch_level_start(1) -        &
360                                            dz(1)/2.0) / dz(1) )               &
361                                      * dz(1) + dz(1)/2.0
362       ENDIF
363       
364       IF ( number_stretch_level_start > 1 ) THEN
365          DO n = 2, number_stretch_level_start
366             dz_stretch_level_start(n) = INT( dz_stretch_level_start(n) /      &
367                                              dz(n) ) * dz(n)
368          ENDDO
369       ENDIF
370       
371       IF ( number_stretch_level_end /= 0 ) THEN
372          DO n = 1, number_stretch_level_end
373             dz_stretch_level_end(n) = INT( dz_stretch_level_end(n) /          &
374                                            dz(n+1) ) * dz(n+1)
375          ENDDO
376       ENDIF
377
378!
379!--    Determine stretching factor if necessary
380       IF ( number_stretch_level_end >= 1 ) THEN
381          CALL calculate_stretching_factor( number_stretch_level_end )
382       ENDIF
383
384!
385!--    Grid for atmosphere with surface at z=0 (k=0, w-grid).
386!--    First compute the u- and v-levels. In case of dirichlet bc for u and v
387!--    the first u/v- and w-level (k=0) are defined at same height (z=0).
388!--    The second u-level (k=1) corresponds to the top of the
389!--    surface layer. In case of symmetric boundaries (closed channel flow),
390!--    the first grid point is always at z=0.
391       IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2 .OR.                              & 
392            topography == 'closed_channel' ) THEN
393          zu(0) = 0.0_wp
394       ELSE
395          zu(0) = - dz(1) * 0.5_wp
396       ENDIF
397         
398       zu(1) =   dz(1) * 0.5_wp
399       
400!
401!--    Determine u and v height levels considering the possibility of grid
402!--    stretching in several heights.
403       n = 1
404       dz_stretch_level_start_index = nzt+1
405       dz_stretch_level_end_index = nzt+1
406       dz_stretched = dz(1)
407
408!--    The default value of dz_stretch_level_start is negative, thus the first
409!--    condition is true even if no stretching shall be applied. Hence, the
410!--    second condition is also necessary.
411       DO  k = 2, nzt+1-symmetry_flag
412          IF ( dz_stretch_level_start(n) <= zu(k-1) .AND.                      &
413               dz_stretch_level_start(n) /= -9999999.9_wp ) THEN
414             dz_stretched = dz_stretched * dz_stretch_factor_array(n)
415             
416             IF ( dz(n) > dz(n+1) ) THEN
417                dz_stretched = MAX( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (higher) dz
418             ELSE
419                dz_stretched = MIN( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (lower) dz
420             ENDIF
421             
422             IF ( dz_stretch_level_start_index(n) == nzt+1 )                         &
423             dz_stretch_level_start_index(n) = k-1
424             
425          ENDIF
426         
427          zu(k) = zu(k-1) + dz_stretched
428         
429!
430!--       Make sure that the stretching ends exactly at dz_stretch_level_end
431          dz_level_end = ABS( zu(k) - dz_stretch_level_end(n) ) 
432         
433          IF ( dz_level_end  < dz(n+1)/3.0 ) THEN
434             zu(k) = dz_stretch_level_end(n)
435             dz_stretched = dz(n+1)
436             dz_stretch_level_end_index(n) = k
437             n = n + 1             
438          ENDIF
439       ENDDO
440       
441!
442!--    If a closed channel flow is simulated, make sure that grid structure is 
443!--    the same for both bottom and top boundary. (Hint: Using a different dz
444!--    at the bottom and at the top makes no sense due to symmetric boundaries
445!--    where dz should be equal. Therefore, different dz at the bottom and top 
446!--    causes an abort (see check_parameters).)
447       IF ( topography == 'closed_channel' ) THEN
448          zu(nzt+1) = zu(nzt) + dz(1) * 0.5_wp
449       ENDIF
450
451!
452!--    Compute the w-levels. They are always staggered half-way between the
453!--    corresponding u-levels. In case of dirichlet bc for u and v at the
454!--    ground the first u- and w-level (k=0) are defined at same height (z=0).
455!--    Per default, the top w-level is extrapolated linearly. In case of
456!--    a closed channel flow, zu(nzt+1) and zw(nzt) must be set explicitely.
457!--    (Hint: Using a different dz at the bottom and at the top makes no sense
458!--    due to symmetric boundaries where dz should be equal. Therefore,
459!--    different dz at the bottom and top causes an abort (see
460!--    check_parameters).)
461       zw(0) = 0.0_wp
462       DO  k = 1, nzt-symmetry_flag
463          zw(k) = ( zu(k) + zu(k+1) ) * 0.5_wp
464       ENDDO
465       IF ( topography == 'closed_channel' ) THEN
466          zw(nzt)   = zw(nzt-1) + dz(1)
467          zw(nzt+1) = zw(nzt) + dz(1)
468       ELSE
469          zw(nzt+1) = zw(nzt) + 2.0_wp * ( zu(nzt+1) - zw(nzt) )
470       ENDIF
471
472    ELSE !ocean branch
473
474!
475!--    The stretching region has to be large enough to allow for a smooth
476!--    transition between two different grid spacings. The number 4 is an
477!--    empirical value
478       DO n = 1, number_stretch_level_start
479          min_dz_stretch_level_end(n) = dz_stretch_level_start(n) -            &
480                                        4 * MAX( dz(n),dz(n+1) )
481       ENDDO
482       
483       IF ( ANY( min_dz_stretch_level_end (1:number_stretch_level_start) <     &
484                 dz_stretch_level_end(1:number_stretch_level_start) ) ) THEN
485             message_string= 'Each dz_stretch_level_end has to be less ' //   &
486                             'than its corresponding value for &' //           &
487                             'dz_stretch_level_start - 4*MAX(dz(n),dz(n+1)) '//&
488                             'to allow for smooth grid stretching'
489             CALL message( 'init_grid', 'PA0224', 1, 2, 0, 6, 0 )
490       ENDIF
491       
492!
493!--    Stretching must not be applied close to the surface (last two grid
494!--    points). For the default case dz_stretch_level_start is negative.
495       IF ( ANY( dz_stretch_level_start >= - dz(1) * 1.5_wp ) ) THEN
496          WRITE( message_string, * ) 'Each dz_stretch_level_start has to be ',&
497                                     'less than ', -dz(1) * 1.5
498             CALL message( 'init_grid', 'PA0226', 1, 2, 0, 6, 0 )
499       ENDIF
500
501!
502!--    The stretching has to start and end on a grid level. Therefore
503!--    user-specified values are mapped to the next highest level. The 
504!--    calculation of the first level is realized differently just because of
505!--    historical reasons (the advanced/new stretching mechanism was realized 
506!--    in a way that results don't change if the old parameters
507!--    dz_stretch_level, dz_stretch_factor and dz_max are used)
508       IF ( number_stretch_level_start /= 0 ) THEN
509          dz_stretch_level_start(1) = INT( (dz_stretch_level_start(1) +        &
510                                            dz(1)/2.0) / dz(1) )               &
511                                      * dz(1) - dz(1)/2.0
512       ENDIF
513       
514       IF ( number_stretch_level_start > 1 ) THEN
515          DO n = 2, number_stretch_level_start
516             dz_stretch_level_start(n) = INT( dz_stretch_level_start(n) /      &
517                                              dz(n) ) * dz(n)
518          ENDDO
519       ENDIF
520       
521       IF ( number_stretch_level_end /= 0 ) THEN
522          DO n = 1, number_stretch_level_end
523             dz_stretch_level_end(n) = INT( dz_stretch_level_end(n) /          &
524                                            dz(n+1) ) * dz(n+1)
525          ENDDO
526       ENDIF
527       
528!
529!--    Determine stretching factor if necessary
530       IF ( number_stretch_level_end >= 1 ) THEN
531          CALL calculate_stretching_factor( number_stretch_level_end )
532       ENDIF
533
534!
535!--    Grid for ocean with free water surface is at k=nzt (w-grid).
536!--    In case of neumann bc at the ground the first first u-level (k=0) lies
537!--    below the first w-level (k=0). In case of dirichlet bc the first u- and
538!--    w-level are defined at same height, but staggered from the second level.
539!--    The second u-level (k=1) corresponds to the top of the surface layer.
540!--    z values are negative starting from z=0 (surface)
541       zu(nzt+1) =   dz(1) * 0.5_wp
542       zu(nzt)   = - dz(1) * 0.5_wp
543
544!
545!--    Determine u and v height levels considering the possibility of grid
546!--    stretching in several heights.
547       n = 1
548       dz_stretch_level_start_index = 0
549       dz_stretch_level_end_index = 0
550       dz_stretched = dz(1)
551
552       DO  k = nzt-1, 0, -1
553         
554          IF ( dz_stretch_level_start(n) >= zu(k+1) ) THEN
555             dz_stretched = dz_stretched * dz_stretch_factor_array(n)
556
557             IF ( dz(n) > dz(n+1) ) THEN
558                dz_stretched = MAX( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (higher) dz
559             ELSE
560                dz_stretched = MIN( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (lower) dz
561             ENDIF
562             
563             IF ( dz_stretch_level_start_index(n) == 0 )                             &
564             dz_stretch_level_start_index(n) = k+1
565             
566          ENDIF
567         
568          zu(k) = zu(k+1) - dz_stretched
569         
570!
571!--       Make sure that the stretching ends exactly at dz_stretch_level_end
572          dz_level_end = ABS( zu(k) - dz_stretch_level_end(n) ) 
573         
574          IF ( dz_level_end  < dz(n+1)/3.0 ) THEN
575             zu(k) = dz_stretch_level_end(n)
576             dz_stretched = dz(n+1)
577             dz_stretch_level_end_index(n) = k
578             n = n + 1             
579          ENDIF
580       ENDDO
581       
582!
583!--    Compute the w-levels. They are always staggered half-way between the
584!--    corresponding u-levels, except in case of dirichlet bc for u and v
585!--    at the ground. In this case the first u- and w-level are defined at
586!--    same height. The top w-level (nzt+1) is not used but set for
587!--    consistency, since w and all scalar variables are defined up tp nzt+1.
588       zw(nzt+1) = dz(1)
589       zw(nzt)   = 0.0_wp
590       DO  k = 0, nzt
591          zw(k) = ( zu(k) + zu(k+1) ) * 0.5_wp
592       ENDDO
593
594!
595!--    In case of dirichlet bc for u and v the first u- and w-level are defined
596!--    at same height.
597       IF ( ibc_uv_b == 0 ) THEN
598          zu(0) = zw(0)
599       ENDIF
600
601    ENDIF !End of defining the vertical grid levels
602
603!
604!-- Compute grid lengths.
605    DO  k = 1, nzt+1
606       dzu(k)  = zu(k) - zu(k-1)
607       ddzu(k) = 1.0_wp / dzu(k)
608       dzw(k)  = zw(k) - zw(k-1)
609       ddzw(k) = 1.0_wp / dzw(k)
610    ENDDO
611
612    DO  k = 1, nzt
613       dd2zu(k) = 1.0_wp / ( dzu(k) + dzu(k+1) )
614    ENDDO
615   
616!   
617!-- The FFT- SOR-pressure solvers assume grid spacings of a staggered grid
618!-- everywhere. For the actual grid, the grid spacing at the lowest level
619!-- is only dz/2, but should be dz. Therefore, an additional array
620!-- containing with appropriate grid information is created for these
621!-- solvers.
622    IF ( psolver(1:9) /= 'multigrid' )  THEN
623       ALLOCATE( ddzu_pres(1:nzt+1) )
624       ddzu_pres = ddzu
625       ddzu_pres(1) = ddzu_pres(2)  ! change for lowest level
626    ENDIF
627
628!
629!-- Compute the reciprocal values of the horizontal grid lengths.
630    ddx = 1.0_wp / dx
631    ddy = 1.0_wp / dy
632    dx2 = dx * dx
633    dy2 = dy * dy
634    ddx2 = 1.0_wp / dx2
635    ddy2 = 1.0_wp / dy2
636
637!
638!-- Allocate 3D array to set topography
639    ALLOCATE( topo(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
640    topo = 0
641!
642!-- Initialize topography by generic topography or read topography from file. 
643    CALL init_topo( topo )
644!
645!-- Set flags to mask topography on the grid.
646    CALL set_topo_flags( topo )
647
648!
649!-- Determine the maximum level of topography. It is used for
650!-- steering the degradation of order of the applied advection scheme,
651!-- as well in the lpm.
652    k_top = 0
653    DO  i = nxl, nxr
654       DO  j = nys, nyn
655          DO  k = nzb, nzt + 1
656             k_top = MAX( k_top, MERGE( k, 0, .NOT. BTEST( topo(k,j,i), 0 ) ) )
657          ENDDO
658       ENDDO
659    ENDDO
660#if defined( __parallel )
661    CALL MPI_ALLREDUCE( k_top, nzb_max, 1, MPI_INTEGER,                        &
662                        MPI_MAX, comm2d, ierr )
663#else
664    nzb_max = k_top
665#endif
666!
667!-- Increment nzb_max by 1 in order to allow for proper diverengence correction.
668!-- Further, in case topography extents up to the model top, limit to nzt.
669    nzb_max = MIN( nzb_max+1, nzt ) 
670!
671!-- Determine minimum index of topography. Usually, this will be nzb. In case
672!-- there is elevated topography, however, the lowest topography will be higher.
673!-- This index is e.g. used to calculate mean first-grid point atmosphere
674!-- temperature, surface pressure and density, etc. .
675    topo_min_level   = 0
676#if defined( __parallel )
677    CALL MPI_ALLREDUCE( MINVAL( topo_top_ind(nys:nyn,nxl:nxr,0) ),             &
678                        topo_min_level, 1, MPI_INTEGER, MPI_MIN, comm2d, ierr )
679#else
680    topo_min_level = MINVAL( topo_top_ind(nys:nyn,nxl:nxr,0) )
681#endif
682
683!
684!-- Check topography for consistency with model domain. Therefore, use
685!-- maximum and minium topography-top indices. Note, minimum topography top
686!-- index is already calculated. 
687    IF ( TRIM( topography ) /= 'flat' )  THEN
688#if defined( __parallel )
689       CALL MPI_ALLREDUCE( MAXVAL( topo_top_ind(nys:nyn,nxl:nxr,0) ),          &
690                           nzb_local_max, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr )               
691#else
692       nzb_local_max = MAXVAL( topo_top_ind(nys:nyn,nxl:nxr,0) )
693#endif
694       nzb_local_min = topo_min_level
695!
696!--    Consistency checks
697       IF ( nzb_local_min < 0  .OR.  nzb_local_max  > nz + 1 )  THEN
698          WRITE( message_string, * ) 'nzb_local values are outside the',       &
699                                ' model domain',                               &
700                                '&MINVAL( nzb_local ) = ', nzb_local_min,      &
701                                '&MAXVAL( nzb_local ) = ', nzb_local_max
702          CALL message( 'init_grid', 'PA0210', 1, 2, 0, 6, 0 )
703       ENDIF
704    ENDIF
705!
706!-- Define vertical gridpoint from (or to) which on the usual finite difference
707!-- form (which does not use surface fluxes) is applied
708    IF ( constant_flux_layer  .OR.  use_surface_fluxes )  THEN
709       nzb_diff = nzb + 2
710    ELSE
711       nzb_diff = nzb + 1
712    ENDIF
713
714    IF ( TRIM( topography ) /= 'flat' )  THEN
715!
716!--    Allocate and set the arrays containing the topography height (for output
717!--    reasons only).
718       IF ( nxr == nx  .AND.  nyn /= ny )  THEN
719          ALLOCATE( zu_s_inner(nxl:nxr+1,nys:nyn),                             &
720                    zw_w_inner(nxl:nxr+1,nys:nyn) )
721       ELSEIF ( nxr /= nx  .AND.  nyn == ny )  THEN
722          ALLOCATE( zu_s_inner(nxl:nxr,nys:nyn+1),                             &
723                    zw_w_inner(nxl:nxr,nys:nyn+1) )
724       ELSEIF ( nxr == nx  .AND.  nyn == ny )  THEN
725          ALLOCATE( zu_s_inner(nxl:nxr+1,nys:nyn+1),                           &
726                    zw_w_inner(nxl:nxr+1,nys:nyn+1) )
727       ELSE
728          ALLOCATE( zu_s_inner(nxl:nxr,nys:nyn),                               &
729                    zw_w_inner(nxl:nxr,nys:nyn) )
730       ENDIF
731
732       zu_s_inner   = 0.0_wp
733       zw_w_inner   = 0.0_wp
734!
735!--    Determine local topography height on scalar and w-grid. Note, setting
736!--    lateral boundary values is not necessary, realized via wall_flags_static_0
737!--    array. Further, please note that loop bounds are different from
738!--    nxl to nxr and nys to nyn on south and right model boundary, hence,
739!--    use intrinsic lbound and ubound functions to infer array bounds.
740       DO  i = LBOUND(zu_s_inner, 1), UBOUND(zu_s_inner, 1)
741          DO  j = LBOUND(zu_s_inner, 2), UBOUND(zu_s_inner, 2)
742!
743!--          Topography height on scalar grid. Therefore, determine index of
744!--          upward-facing surface element on scalar grid.
745             zu_s_inner(i,j) = zu(topo_top_ind(j,i,0))
746!
747!--          Topography height on w grid. Therefore, determine index of
748!--          upward-facing surface element on w grid.
749             zw_w_inner(i,j) = zw(topo_top_ind(j,i,3))
750          ENDDO
751       ENDDO
752    ENDIF
753!
754!-- Vertical nesting: communicate vertical grid level arrays between fine and
755!-- coarse grid
756    IF ( vnested )  CALL vnest_init_grid
757
758 END SUBROUTINE init_grid
759
760
761! Description:
762! -----------------------------------------------------------------------------!
763!> Calculation of the stretching factor through an iterative method. Ideas were
764!> taken from the paper "Regional stretched grid generation and its application
765!> to the NCAR RegCM (1999)". Normally, no analytic solution exists because the
766!> system of equations has two variables (r,l) but four requirements
767!> (l=integer, r=[0,88;1,2], Eq(6), Eq(5) starting from index j=1) which
768!> results into an overdetermined system.
769!------------------------------------------------------------------------------!
770 SUBROUTINE calculate_stretching_factor( number_end )
771 
772    USE control_parameters,                                                    &
773        ONLY:  dz, dz_stretch_factor_array,                 &
774               dz_stretch_level_end, dz_stretch_level_start, message_string
775 
776    USE kinds
777   
778    IMPLICIT NONE
779   
780    INTEGER(iwp) ::  iterations  !< number of iterations until stretch_factor_lower/upper_limit is reached 
781    INTEGER(iwp) ::  l_rounded   !< after l_rounded grid levels dz(n) is strechted to dz(n+1) with stretch_factor_2
782    INTEGER(iwp) ::  n           !< loop variable for stretching
783   
784    INTEGER(iwp), INTENT(IN) ::  number_end !< number of user-specified end levels for stretching
785       
786    REAL(wp) ::  delta_l               !< absolute difference between l and l_rounded
787    REAL(wp) ::  delta_stretch_factor  !< absolute difference between stretch_factor_1 and stretch_factor_2
788    REAL(wp) ::  delta_total_new       !< sum of delta_l and delta_stretch_factor for the next iteration (should be as small as possible)
789    REAL(wp) ::  delta_total_old       !< sum of delta_l and delta_stretch_factor for the last iteration
790    REAL(wp) ::  distance              !< distance between dz_stretch_level_start and dz_stretch_level_end (stretching region)
791    REAL(wp) ::  l                     !< value that fulfil Eq. (5) in the paper mentioned above together with stretch_factor_1 exactly
792    REAL(wp) ::  numerator             !< numerator of the quotient
793    REAL(wp) ::  stretch_factor_1      !< stretching factor that fulfil Eq. (5) togehter with l exactly
794    REAL(wp) ::  stretch_factor_2      !< stretching factor that fulfil Eq. (6) togehter with l_rounded exactly
795   
796    REAL(wp) ::  dz_stretch_factor_array_2(9) = 1.08_wp  !< Array that contains all stretch_factor_2 that belongs to stretch_factor_1
797   
798    REAL(wp), PARAMETER ::  stretch_factor_interval = 1.0E-06_wp  !< interval for sampling possible stretching factors
799    REAL(wp), PARAMETER ::  stretch_factor_lower_limit = 0.88_wp  !< lowest possible stretching factor
800    REAL(wp), PARAMETER ::  stretch_factor_upper_limit = 1.12_wp  !< highest possible stretching factor
801 
802 
803    l = 0
804    DO  n = 1, number_end
805   
806       iterations = 1
807       stretch_factor_1 = 1.0_wp 
808       stretch_factor_2 = 1.0_wp
809       delta_total_old = 1.0_wp
810       
811!
812!--    First branch for stretching from rough to fine
813       IF ( dz(n) > dz(n+1) ) THEN
814          DO WHILE ( stretch_factor_1 >= stretch_factor_lower_limit ) 
815             
816             stretch_factor_1 = 1.0_wp - iterations * stretch_factor_interval
817             distance = ABS( dz_stretch_level_end(n) -                         &
818                        dz_stretch_level_start(n) )   
819             numerator = distance*stretch_factor_1/dz(n) +                     &
820                         stretch_factor_1 - distance/dz(n)
821             
822             IF ( numerator > 0.0_wp ) THEN
823                l = LOG( numerator ) / LOG( stretch_factor_1 ) - 1.0_wp
824                l_rounded = NINT( l )
825                delta_l = ABS( l_rounded - l ) / l
826             ENDIF
827             
828             stretch_factor_2 = EXP( LOG( dz(n+1)/dz(n) ) / (l_rounded) )
829             
830             delta_stretch_factor = ABS( stretch_factor_1 -                    &
831                                         stretch_factor_2 ) /                  &
832                                    stretch_factor_2
833             
834             delta_total_new = delta_l + delta_stretch_factor
835
836!
837!--          stretch_factor_1 is taken to guarantee that the stretching
838!--          procedure ends as close as possible to dz_stretch_level_end.
839!--          stretch_factor_2 would guarantee that the stretched dz(n) is
840!--          equal to dz(n+1) after l_rounded grid levels.
841             IF (delta_total_new < delta_total_old) THEN
842                dz_stretch_factor_array(n) = stretch_factor_1
843                dz_stretch_factor_array_2(n) = stretch_factor_2
844                delta_total_old = delta_total_new
845             ENDIF
846             
847             iterations = iterations + 1
848           
849          ENDDO
850
851!
852!--    Second branch for stretching from fine to rough
853       ELSEIF ( dz(n) < dz(n+1) ) THEN
854          DO WHILE ( stretch_factor_1 <= stretch_factor_upper_limit )
855                     
856             stretch_factor_1 = 1.0_wp + iterations * stretch_factor_interval
857             distance = ABS( dz_stretch_level_end(n) -                         &
858                        dz_stretch_level_start(n) ) 
859             numerator = distance*stretch_factor_1/dz(n) +                     &
860                         stretch_factor_1 - distance/dz(n)
861             
862             l = LOG( numerator ) / LOG( stretch_factor_1 ) - 1.0_wp
863             l_rounded = NINT( l )
864             delta_l = ABS( l_rounded - l ) / l
865             
866             stretch_factor_2 = EXP( LOG( dz(n+1)/dz(n) ) / (l_rounded) )
867
868             delta_stretch_factor = ABS( stretch_factor_1 -                    &
869                                        stretch_factor_2 ) /                   &
870                                        stretch_factor_2
871             
872             delta_total_new = delta_l + delta_stretch_factor
873             
874!
875!--          stretch_factor_1 is taken to guarantee that the stretching
876!--          procedure ends as close as possible to dz_stretch_level_end.
877!--          stretch_factor_2 would guarantee that the stretched dz(n) is
878!--          equal to dz(n+1) after l_rounded grid levels.
879             IF (delta_total_new < delta_total_old) THEN
880                dz_stretch_factor_array(n) = stretch_factor_1
881                dz_stretch_factor_array_2(n) = stretch_factor_2
882                delta_total_old = delta_total_new
883             ENDIF
884             
885             iterations = iterations + 1
886          ENDDO
887         
888       ELSE
889          message_string= 'Two adjacent values of dz must be different'
890          CALL message( 'init_grid', 'PA0228', 1, 2, 0, 6, 0 )
891         
892       ENDIF
893
894!
895!--    Check if also the second stretching factor fits into the allowed
896!--    interval. If not, print a warning for the user.
897       IF ( dz_stretch_factor_array_2(n) < stretch_factor_lower_limit .OR.     & 
898            dz_stretch_factor_array_2(n) > stretch_factor_upper_limit ) THEN
899          WRITE( message_string, * ) 'stretch_factor_2 = ',                    &
900                                     dz_stretch_factor_array_2(n), ' which is',&
901                                     ' responsible for exactly reaching& dz =',&
902                                      dz(n+1), 'after a specific amount of',   & 
903                                     ' grid levels& exceeds the upper',        &
904                                     ' limit =', stretch_factor_upper_limit,   &
905                                     ' &or lower limit = ',                    &
906                                     stretch_factor_lower_limit
907          CALL message( 'init_grid', 'PA0499', 0, 1, 0, 6, 0 )
908           
909       ENDIF
910    ENDDO
911       
912 END SUBROUTINE calculate_stretching_factor
913 
914 
915! Description:
916! -----------------------------------------------------------------------------!
917!> Set temporary topography flags and reference buildings on top of underlying
918!> orography.
919!------------------------------------------------------------------------------!
920 SUBROUTINE process_topography( topo_3d )
921
922    USE arrays_3d,                                                             &
923        ONLY:  zu, zw
924
925    USE control_parameters,                                                    &
926        ONLY:  bc_lr_cyc, bc_ns_cyc, message_string, ocean_mode
927
928    USE indices,                                                               &
929        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb,  &
930               nzt
931
932    USE netcdf_data_input_mod,                                                 &
933        ONLY:  buildings_f, building_id_f, building_type_f, input_pids_static, &
934               terrain_height_f
935
936    USE kinds
937
938    USE pegrid
939
940    IMPLICIT NONE
941
942    INTEGER(iwp) ::  i                !< running index along x-direction
943    INTEGER(iwp) ::  j                !< running index along y-direction
944    INTEGER(iwp) ::  k                !< running index along z-direction with respect to numeric grid
945    INTEGER(iwp) ::  k2               !< running index along z-direction with respect to netcdf grid
946    INTEGER(iwp) ::  nr               !< index variable indication maximum terrain height for respective building ID
947    INTEGER(iwp) ::  num_build        !< counter for number of buildings
948    INTEGER(iwp) ::  topo_top_index   !< orography top index, used to map 3D buildings onto terrain
949
950    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  displace_dum        !< displacements of start addresses, used for MPI_ALLGATHERV
951    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  build_ids           !< building IDs on entire model domain
952    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  build_ids_final     !< building IDs on entire model domain, multiple occurences are sorted out
953    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  build_ids_final_tmp !< temporary array used for resizing
954    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  build_ids_l         !< building IDs on local subdomain
955    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  build_ids_l_tmp     !< temporary array used to resize array of building IDs
956
957    INTEGER(iwp), DIMENSION(0:numprocs-1) ::  num_buildings     !< number of buildings with different ID on entire model domain
958    INTEGER(iwp), DIMENSION(0:numprocs-1) ::  num_buildings_l   !< number of buildings with different ID on local subdomain
959
960    INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  topo_3d !< input array for 3D topography and dummy array for setting "outer"-flags
961
962    REAL(wp)                            ::  ocean_offset        !< offset to consider inverse vertical coordinate at topography definition
963    REAL(wp)                            ::  oro_min = 0.0_wp    !< minimum terrain height in entire model domain, used to reference terrain to zero
964    REAL(wp), DIMENSION(:), ALLOCATABLE ::  oro_max             !< maximum terrain height occupied by an building with certain id
965    REAL(wp), DIMENSION(:), ALLOCATABLE ::  oro_max_l           !< maximum terrain height occupied by an building with certain id, on local subdomain
966
967!
968!-- Reference lowest terrain height to zero. In case the minimum terrain height
969!-- is non-zero, all grid points of the lower vertical grid levels might be
970!-- entirely below the surface, meaning a waste of computational resources.
971!-- In order to avoid this, remove the lowest terrain height. Please note,
972!-- in case of a nested run, the global minimum from all parent and childs
973!-- need to be remove to avoid steep edges at the child-domain boundaries.
974    IF ( input_pids_static )  THEN
975   
976#if defined( __parallel ) 
977       CALL MPI_ALLREDUCE( MINVAL( terrain_height_f%var ), oro_min, 1,         &
978                           MPI_REAL, MPI_MIN, MPI_COMM_WORLD, ierr )
979#else
980       oro_min = MINVAL( terrain_height_f%var )
981#endif
982       terrain_height_f%var = terrain_height_f%var - oro_min
983!                           
984!--    Give an informative message that terrain height is referenced to zero   
985       IF ( oro_min > 0.0_wp )  THEN
986          WRITE( message_string, * ) 'Terrain height was internally shifted '//&
987                          'downwards by ', oro_min, 'meter(s) to save ' //     &
988                          'computational resources.'
989          CALL message( 'init_grid', 'PA0505', 0, 0, 0, 6, 0 )
990       ENDIF
991    ENDIF   
992   
993!
994!-- In the following, buildings and orography are further preprocessed
995!-- before they are mapped on the LES grid.
996!-- Buildings are mapped on top of the orography by maintaining the roof
997!-- shape of the building. This can be achieved by referencing building on
998!-- top of the maximum terrain height within the area occupied by the
999!-- respective building. As buildings and terrain height are defined PE-wise,
1000!-- parallelization of this referencing is required (a building can be
1001!-- distributed between different PEs). 
1002!-- In a first step, determine the number of buildings with different
1003!-- building id on each PE. In a next step, all building ids are gathered
1004!-- into one array which is present to all PEs. For each building ID,
1005!-- the maximum terrain height occupied by the respective building is
1006!-- computed and distributed to each PE. 
1007!-- Finally, for each building id and its respective reference orography,
1008!-- builidings are mapped on top.   
1009!--
1010!-- First, pre-set topography flags, bit 1 indicates orography, bit 2
1011!-- buildings
1012!-- classify the respective surfaces.
1013    topo_3d          = IBSET( topo_3d, 0 )
1014    topo_3d(nzb,:,:) = IBCLR( topo_3d(nzb,:,:), 0 )
1015!
1016!-- In order to map topography on PALM grid also in case of ocean simulations,
1017!-- pre-calculate an offset value.
1018    ocean_offset = MERGE( zw(0), 0.0_wp, ocean_mode )
1019!
1020!-- Reference buildings on top of orography. This is not necessary
1021!-- if topography is read from ASCII file as no distinction between buildings
1022!-- and terrain height can be made. Moreover, this is also not necessary if
1023!-- urban-surface and land-surface model are used at the same time.
1024    IF ( input_pids_static )  THEN
1025
1026       IF ( buildings_f%from_file )  THEN
1027          num_buildings_l = 0
1028          num_buildings   = 0
1029!
1030!--       Allocate at least one element for building ids and give it an inital
1031!--       negative value that will be overwritten later. This, however, is
1032!--       necessary in case there all IDs in the model domain are fill values.
1033          ALLOCATE( build_ids_l(1) )
1034          build_ids_l = -1 
1035          DO  i = nxl, nxr
1036             DO  j = nys, nyn
1037                IF ( building_id_f%var(j,i) /= building_id_f%fill )  THEN
1038                   IF ( num_buildings_l(myid) > 0 )  THEN
1039                      IF ( ANY( building_id_f%var(j,i) ==  build_ids_l ) )   &
1040                      THEN
1041                         CYCLE
1042                      ELSE
1043                         num_buildings_l(myid) = num_buildings_l(myid) + 1
1044!
1045!--                   Resize array with different local building ids
1046                      ALLOCATE( build_ids_l_tmp(1:SIZE(build_ids_l)) )
1047                      build_ids_l_tmp = build_ids_l
1048                      DEALLOCATE( build_ids_l )
1049                      ALLOCATE( build_ids_l(1:num_buildings_l(myid)) )
1050                      build_ids_l(1:num_buildings_l(myid)-1) =                 &
1051                                  build_ids_l_tmp(1:num_buildings_l(myid)-1)
1052                      build_ids_l(num_buildings_l(myid)) = building_id_f%var(j,i)
1053                      DEALLOCATE( build_ids_l_tmp )
1054                   ENDIF
1055!
1056!--                First occuring building id on PE
1057                   ELSE
1058                      num_buildings_l(myid) = num_buildings_l(myid) + 1
1059                      build_ids_l(1) = building_id_f%var(j,i)
1060                   ENDIF
1061                ENDIF
1062             ENDDO
1063          ENDDO
1064!
1065!--       Determine number of different building ids for the entire domain
1066#if defined( __parallel ) 
1067          CALL MPI_ALLREDUCE( num_buildings_l, num_buildings, numprocs,              &
1068                              MPI_INTEGER, MPI_SUM, comm2d, ierr ) 
1069#else
1070          num_buildings = num_buildings_l
1071#endif
1072!
1073!--       Gather all buildings ids on each PEs.
1074!--       First, allocate array encompassing all building ids in model domain. 
1075          ALLOCATE( build_ids(1:SUM(num_buildings)) )
1076#if defined( __parallel ) 
1077!
1078!--       Allocate array for displacements.
1079!--       As each PE may has a different number of buildings, so that
1080!--       the block sizes send by each PE may not be equal. Hence,
1081!--       information about the respective displacement is required, indicating
1082!--       the respective adress where each MPI-task writes into the receive
1083!--       buffer array 
1084          ALLOCATE( displace_dum(0:numprocs-1) )
1085          displace_dum(0) = 0
1086          DO i = 1, numprocs-1
1087             displace_dum(i) = displace_dum(i-1) + num_buildings(i-1)
1088          ENDDO
1089
1090          CALL MPI_ALLGATHERV( build_ids_l(1:num_buildings_l(myid)),                 &
1091                               num_buildings(myid),                                  &
1092                               MPI_INTEGER,                                          &
1093                               build_ids,                                            &
1094                               num_buildings,                                        &
1095                               displace_dum,                                         & 
1096                               MPI_INTEGER,                                          &
1097                               comm2d, ierr )   
1098
1099          DEALLOCATE( displace_dum )
1100
1101#else
1102          build_ids = build_ids_l
1103#endif
1104
1105!
1106!--       Note, in parallel mode building ids can occure mutliple times, as
1107!--       each PE has send its own ids. Therefore, sort out building ids which
1108!--       appear more than one time.
1109          num_build = 0
1110          DO  nr = 1, SIZE(build_ids)
1111
1112             IF ( ALLOCATED(build_ids_final) )  THEN
1113                IF ( ANY( build_ids(nr) == build_ids_final ) )  THEN
1114                   CYCLE
1115                ELSE
1116                   num_build = num_build + 1
1117!
1118!--                Resize
1119                   ALLOCATE( build_ids_final_tmp(1:num_build) )
1120                   build_ids_final_tmp(1:num_build-1) = build_ids_final(1:num_build-1)
1121                   DEALLOCATE( build_ids_final )
1122                   ALLOCATE( build_ids_final(1:num_build) )
1123                   build_ids_final(1:num_build-1) = build_ids_final_tmp(1:num_build-1)
1124                   build_ids_final(num_build) = build_ids(nr)
1125                   DEALLOCATE( build_ids_final_tmp )
1126                ENDIF             
1127             ELSE
1128                num_build = num_build + 1
1129                ALLOCATE( build_ids_final(1:num_build) )
1130                build_ids_final(num_build) = build_ids(nr)
1131             ENDIF
1132          ENDDO
1133
1134!
1135!--       Determine maximumum terrain height occupied by the respective
1136!--       building and temporalily store on oro_max
1137          ALLOCATE( oro_max_l(1:SIZE(build_ids_final)) )
1138          ALLOCATE( oro_max(1:SIZE(build_ids_final))   )
1139          oro_max_l = 0.0_wp
1140
1141          DO  nr = 1, SIZE(build_ids_final)
1142             oro_max_l(nr) = MAXVAL(                                           &
1143                              MERGE( terrain_height_f%var(nys:nyn,nxl:nxr),    &
1144                                     0.0_wp,                                   &
1145                                     building_id_f%var(nys:nyn,nxl:nxr) ==     &
1146                                     build_ids_final(nr) ) )
1147          ENDDO
1148   
1149#if defined( __parallel )   
1150          IF ( SIZE(build_ids_final) >= 1 ) THEN
1151             CALL MPI_ALLREDUCE( oro_max_l, oro_max, SIZE( oro_max ), MPI_REAL,&
1152                                 MPI_MAX, comm2d, ierr ) 
1153          ENDIF
1154#else
1155          oro_max = oro_max_l
1156#endif
1157!
1158!--       Finally, determine discrete grid height of maximum orography occupied
1159!--       by a building. Use all-or-nothing approach, i.e. if terrain
1160!--       exceeds the scalar level the grid box is fully terrain and the
1161!--       maximum terrain is set to the zw level.
1162!--       terrain or
1163          oro_max_l = 0.0
1164          DO  nr = 1, SIZE(build_ids_final)
1165             DO  k = nzb, nzt
1166                IF ( zu(k) - ocean_offset <= oro_max(nr) )                     &
1167                   oro_max_l(nr) = zw(k) - ocean_offset
1168             ENDDO
1169             oro_max(nr) = oro_max_l(nr)
1170          ENDDO
1171       ENDIF
1172!
1173!--    Allocate array for storing terrain height under buildings
1174       IF ( buildings_f%from_file )  THEN
1175          ALLOCATE( buildings_f%oro_max(nysg:nyng,nxlg:nxrg) )
1176          buildings_f%oro_max = buildings_f%fill1
1177       END IF
1178!
1179!--    Map orography as well as buildings onto grid.
1180       DO  i = nxl, nxr
1181          DO  j = nys, nyn
1182             topo_top_index = 0
1183!
1184!--          Obtain index in global building_id array
1185             IF ( buildings_f%from_file )  THEN
1186                IF ( building_id_f%var(j,i) /= building_id_f%fill )  THEN
1187!
1188!--                Determine index where maximum terrain height occupied by
1189!--                the respective building height is stored.
1190                   nr = MINLOC( ABS( build_ids_final -                         &
1191                                     building_id_f%var(j,i) ), DIM = 1 )
1192!
1193!--                Save grid-indexed oro_max
1194                   buildings_f%oro_max(j,i) = oro_max(nr)
1195                ENDIF
1196             ENDIF
1197             DO  k = nzb, nzt
1198!
1199!--             In a first step, if grid point is below or equal the given
1200!--             terrain height, grid point is flagged to be of type natural.
1201!--             Please note, in case there is also a building which is lower
1202!--             than the vertical grid spacing, initialization of surface
1203!--             attributes will not be correct as given surface information
1204!--             will not be in accordance to the classified grid points.
1205!--             Hence, in this case, also a building flag.
1206                IF ( zu(k) - ocean_offset <= terrain_height_f%var(j,i) )  THEN
1207                   topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
1208                   topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 1 )
1209                   topo_top_index = k ! topo_top_index + 1
1210                ENDIF
1211!
1212!--             Set building grid points. Here, only consider 2D buildings.
1213!--             3D buildings require separate treatment.
1214                IF ( buildings_f%from_file  .AND.  buildings_f%lod == 1 )  THEN
1215!
1216!--                Fill-up the terrain to the level of maximum orography
1217!--                within the building-covered area.
1218                   IF ( building_id_f%var(j,i) /= building_id_f%fill )  THEN
1219!
1220!--                   Note, oro_max is always on zw level                   
1221                      IF ( zu(k) - ocean_offset < oro_max(nr) )  THEN
1222                         topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
1223                         topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 1 )
1224                      ELSEIF ( zu(k) - ocean_offset <=                         &
1225                               oro_max(nr) + buildings_f%var_2d(j,i) )  THEN
1226                         topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
1227                         topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 2 )
1228                      ENDIF
1229                   ENDIF
1230                ENDIF
1231             ENDDO
1232!
1233!--          Special treatment for non grid-resolved buildings. This case,
1234!--          the uppermost terrain grid point is flagged as building as well
1235!--          well, even though no building exists at all. However, the
1236!--          surface element will be identified as urban-surface and the
1237!--          input data provided by the drivers is consistent to the surface
1238!--          classification. Else, all non grid-resolved buildings would vanish
1239!--          and identified as terrain grid points, which, however, won't be
1240!--          consistent with the input data.
1241             IF ( buildings_f%from_file  .AND.  buildings_f%lod == 1 )  THEN
1242                IF ( building_id_f%var(j,i) /= building_id_f%fill )  THEN
1243                   DO  k = nzb, nzt
1244                      IF( zw(k) - ocean_offset == oro_max(nr) )  THEN
1245                         IF ( buildings_f%var_2d(j,i) <= zu(k+1) - zw(k) )  THEN
1246                            topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 2 )
1247                         ENDIF
1248                      ENDIF
1249                   ENDDO
1250                ENDIF
1251             ENDIF
1252!
1253!--          Map 3D buildings onto terrain height. 
1254!--          In case of any slopes, map building on top of maximum terrain
1255!--          height covered by the building. In other words, extend
1256!--          building down to the respective local terrain-surface height.
1257             IF ( buildings_f%from_file  .AND.  buildings_f%lod == 2 )  THEN
1258                IF ( building_id_f%var(j,i) /= building_id_f%fill )  THEN
1259!
1260!--                Extend building down to the terrain surface, i.e. fill-up
1261!--                surface irregularities below a building. Note, oro_max
1262!--                is already a discrete height according to the all-or-nothing
1263!--                approach, i.e. grid box is either topography or atmosphere,
1264!--                terrain top is defined at upper bound of the grid box.
1265!--                Hence, check for zw in this case.
1266!--                Note, do this only for buildings which are surface mounted,
1267!--                i.e. building types 1-6. Below bridges, which are represented
1268!--                exclusively by building type 7, terrain shape should be
1269!--                maintained.
1270                   IF ( building_type_f%from_file )  THEN
1271                      IF ( building_type_f%var(j,i) /= 7 )  THEN
1272                         DO k = topo_top_index + 1, nzt + 1     
1273                            IF ( zu(k) - ocean_offset <= oro_max(nr) )  THEN
1274                               topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
1275                               topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 1 )
1276                            ENDIF
1277                         ENDDO       
1278!                     
1279!--                      After surface irregularities are smoothen, determine
1280!--                      lower start index where building starts.
1281                         DO  k = nzb, nzt
1282                            IF ( zu(k) - ocean_offset <= oro_max(nr) )         &
1283                               topo_top_index = k
1284                         ENDDO
1285                      ENDIF
1286                   ENDIF
1287!
1288!--                Finally, map building on top.
1289                   k2 = 0
1290                   DO k = topo_top_index, nzt + 1
1291                      IF ( k2 <= buildings_f%nz-1 )  THEN
1292                         IF ( buildings_f%var_3d(k2,j,i) == 1 )  THEN
1293                            topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
1294                            topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 2 )
1295                         ENDIF
1296                      ENDIF
1297                      k2 = k2 + 1
1298                   ENDDO
1299                ENDIF
1300             ENDIF
1301          ENDDO
1302       ENDDO
1303!
1304!--    Horizontal exchange the oro_max array, which is required to for
1305!--    initialization of building-surface properties.
1306       IF ( ALLOCATED( buildings_f%oro_max ) )  THEN
1307          CALL exchange_horiz_2d( buildings_f%oro_max(:,:), nbgp )
1308       ENDIF
1309!
1310!--    Deallocate temporary arrays required for processing and reading data
1311       IF ( ALLOCATED( oro_max         ) )  DEALLOCATE( oro_max         )
1312       IF ( ALLOCATED( oro_max_l       ) )  DEALLOCATE( oro_max_l       )
1313       IF ( ALLOCATED( build_ids_final ) )  DEALLOCATE( build_ids_final )
1314!
1315!-- Topography input via ASCII format.
1316    ELSE
1317       ocean_offset     = MERGE( zw(0), 0.0_wp, ocean_mode )
1318!
1319!--    Initialize topography bit 0 (indicates obstacle) everywhere to zero
1320!--    and clear all grid points at nzb, where alway a surface is defined.
1321!--    Further, set also bit 1 (indicates terrain) at nzb, which is further
1322!--    used for masked data output and further processing. Note, in the
1323!--    ASCII case no distinction is made between buildings and terrain,
1324!--    so that setting of bit 1 and 2 at the same time has no effect.
1325       topo_3d          = IBSET( topo_3d, 0 )
1326       topo_3d(nzb,:,:) = IBCLR( topo_3d(nzb,:,:), 0 )
1327       topo_3d(nzb,:,:) = IBSET( topo_3d(nzb,:,:), 1 )
1328       DO  i = nxl, nxr
1329          DO  j = nys, nyn
1330             DO  k = nzb, nzt
1331!
1332!--             Flag topography for all grid points which are below
1333!--             the local topography height.
1334!--             Note, each topography is flagged as building.
1335                IF ( zu(k) - ocean_offset <= buildings_f%var_2d(j,i) )  THEN
1336                   topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
1337                   topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 2 ) !indicates building
1338                ENDIF
1339             ENDDO
1340          ENDDO
1341       ENDDO
1342    ENDIF
1343
1344    CALL exchange_horiz_int( topo_3d, nys, nyn, nxl, nxr, nzt, nbgp )
1345
1346    IF ( .NOT. bc_ns_cyc )  THEN
1347       IF ( nys == 0  )  topo_3d(:,-1,:)   = topo_3d(:,0,:)
1348       IF ( nyn == ny )  topo_3d(:,ny+1,:) = topo_3d(:,ny,:)
1349    ENDIF
1350
1351    IF ( .NOT. bc_lr_cyc )  THEN
1352       IF ( nxl == 0  )  topo_3d(:,:,-1)   = topo_3d(:,:,0)
1353       IF ( nxr == nx )  topo_3d(:,:,nx+1) = topo_3d(:,:,nx)         
1354    ENDIF
1355
1356 END SUBROUTINE process_topography
1357
1358
1359! Description:
1360! -----------------------------------------------------------------------------!
1361!> Filter topography, i.e. fill holes resolved by only one grid point. 
1362!> Such holes are suspected to lead to velocity blow-ups as continuity
1363!> equation on discrete grid cannot be fulfilled in such case.
1364!------------------------------------------------------------------------------!
1365 SUBROUTINE filter_topography( topo_3d )
1366
1367    USE control_parameters,                                                    &
1368        ONLY:  bc_lr_cyc, bc_ns_cyc, message_string
1369
1370    USE indices,                                                               &
1371        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb, nzt
1372
1373    USE netcdf_data_input_mod,                                                 &
1374        ONLY:  building_id_f, building_type_f 
1375
1376    USE  pegrid
1377
1378    IMPLICIT NONE
1379
1380    LOGICAL      ::  filled = .FALSE. !< flag indicating if holes were filled
1381
1382    INTEGER(iwp) ::  i          !< running index along x-direction
1383    INTEGER(iwp) ::  j          !< running index along y-direction
1384    INTEGER(iwp) ::  k          !< running index along z-direction
1385    INTEGER(iwp) ::  num_hole   !< number of holes (in topography) resolved by only one grid point
1386    INTEGER(iwp) ::  num_hole_l !< number of holes (in topography) resolved by only one grid point on local PE     
1387    INTEGER(iwp) ::  num_wall   !< number of surrounding vertical walls for a single grid point
1388
1389    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE            ::  topo_tmp          !< temporary 3D-topography used to fill holes
1390    INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  topo_3d           !< 3D-topography array merging buildings and orography
1391!
1392!-- Before checking for holes, set lateral boundary conditions for
1393!-- topography. After hole-filling, boundary conditions must be set again.
1394!-- Several iterations are performed, in order to fill holes which might
1395!-- emerge by the filling-algorithm itself.
1396    ALLOCATE( topo_tmp(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1397    topo_tmp = 0
1398
1399    num_hole = 99999
1400    DO WHILE ( num_hole > 0 )       
1401
1402       num_hole = 0   
1403       CALL exchange_horiz_int( topo_3d, nys, nyn, nxl, nxr, nzt, nbgp )
1404!
1405!--    Exchange also building ID and type. Note, building_type is an one-byte
1406!--    variable.
1407       IF ( building_id_f%from_file )                                          &
1408          CALL exchange_horiz_2d_int( building_id_f%var, nys, nyn, nxl, nxr, nbgp )
1409       IF ( building_type_f%from_file )                                        &
1410          CALL exchange_horiz_2d_byte( building_type_f%var, nys, nyn, nxl, nxr, nbgp )
1411
1412       topo_tmp = topo_3d
1413!
1414!--    In case of non-cyclic lateral boundaries, assume lateral boundary to be
1415!--    a solid wall. Thus, intermediate spaces of one grid point between
1416!--    boundary and some topographic structure will be filled.           
1417       IF ( .NOT. bc_ns_cyc )  THEN
1418          IF ( nys == 0  )  topo_tmp(:,-1,:)   = IBCLR( topo_tmp(:,0,:),  0 )
1419          IF ( nyn == ny )  topo_tmp(:,ny+1,:) = IBCLR( topo_tmp(:,ny,:), 0 )
1420       ENDIF
1421
1422       IF ( .NOT. bc_lr_cyc )  THEN
1423          IF ( nxl == 0  )  topo_tmp(:,:,-1)   = IBCLR( topo_tmp(:,:,0),  0 )
1424          IF ( nxr == nx )  topo_tmp(:,:,nx+1) = IBCLR( topo_tmp(:,:,nx), 0 )         
1425       ENDIF
1426
1427       num_hole_l = 0
1428       DO i = nxl, nxr
1429          DO j = nys, nyn
1430             DO  k = nzb+1, nzt
1431                IF ( BTEST( topo_tmp(k,j,i), 0 ) )  THEN
1432                   num_wall = 0
1433                   IF ( .NOT. BTEST( topo_tmp(k,j-1,i), 0 ) )                  &
1434                      num_wall = num_wall + 1
1435                   IF ( .NOT. BTEST( topo_tmp(k,j+1,i), 0 ) )                  &
1436                      num_wall = num_wall + 1
1437                   IF ( .NOT. BTEST( topo_tmp(k,j,i-1), 0 ) )                  &
1438                      num_wall = num_wall + 1
1439                   IF ( .NOT. BTEST( topo_tmp(k,j,i+1), 0 ) )                  &
1440                      num_wall = num_wall + 1
1441                   IF ( .NOT. BTEST( topo_tmp(k-1,j,i), 0 ) )                  &
1442                      num_wall = num_wall + 1   
1443                   IF ( .NOT. BTEST( topo_tmp(k+1,j,i), 0 ) )                  &
1444                      num_wall = num_wall + 1
1445
1446                   IF ( num_wall >= 4 )  THEN
1447                      num_hole_l     = num_hole_l + 1
1448!
1449!--                   Clear flag 0 and set special flag ( bit 4) to indicate
1450!--                   that new topography point is a result of filtering process.
1451                      topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
1452                      topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 4 )
1453!
1454!--                   If filled grid point is occupied by a building, classify
1455!--                   it as building grid point.
1456                      IF ( building_type_f%from_file )  THEN
1457                         IF ( building_type_f%var(j,i)   /=                    & 
1458                              building_type_f%fill            .OR.             &       
1459                              building_type_f%var(j+1,i) /=                    & 
1460                              building_type_f%fill            .OR.             &               
1461                              building_type_f%var(j-1,i) /=                    &               
1462                              building_type_f%fill            .OR.             &               
1463                              building_type_f%var(j,i+1) /=                    &               
1464                              building_type_f%fill            .OR.             &               
1465                              building_type_f%var(j,i-1) /=                    &               
1466                              building_type_f%fill )  THEN
1467!
1468!--                         Set flag indicating building surfaces
1469                            topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 2 )
1470!
1471!--                         Set building_type and ID at this position if not
1472!--                         already set. This is required for proper
1473!--                         initialization of urban-surface energy balance
1474!--                         solver.
1475                            IF ( building_type_f%var(j,i) ==                   &
1476                                 building_type_f%fill )  THEN
1477
1478                               IF ( building_type_f%var(j+1,i) /=              &
1479                                    building_type_f%fill )  THEN
1480                                  building_type_f%var(j,i) =                   &
1481                                                    building_type_f%var(j+1,i)
1482                                  building_id_f%var(j,i) =                     &
1483                                                    building_id_f%var(j+1,i)
1484                               ELSEIF ( building_type_f%var(j-1,i) /=          &
1485                                        building_type_f%fill )  THEN
1486                                  building_type_f%var(j,i) =                   &
1487                                                    building_type_f%var(j-1,i)
1488                                  building_id_f%var(j,i) =                     &
1489                                                    building_id_f%var(j-1,i)
1490                               ELSEIF ( building_type_f%var(j,i+1) /=          &
1491                                        building_type_f%fill )  THEN
1492                                  building_type_f%var(j,i) =                   &
1493                                                    building_type_f%var(j,i+1)
1494                                  building_id_f%var(j,i) =                     &
1495                                                    building_id_f%var(j,i+1)
1496                               ELSEIF ( building_type_f%var(j,i-1) /=          &
1497                                        building_type_f%fill )  THEN
1498                                  building_type_f%var(j,i) =                   &
1499                                                    building_type_f%var(j,i-1)
1500                                  building_id_f%var(j,i) =                     &
1501                                                    building_id_f%var(j,i-1)
1502                               ENDIF
1503                            ENDIF
1504                         ENDIF
1505                      ENDIF
1506!
1507!--                   If filled grid point is already classified as building
1508!--                   everything is fine, else classify this grid point as
1509!--                   natural type grid point. This case, values for the
1510!--                   surface type are already set.
1511                      IF ( .NOT. BTEST( topo_3d(k,j,i), 2 ) )  THEN
1512                         topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 1 )
1513                      ENDIF
1514                   ENDIF
1515                ENDIF
1516             ENDDO
1517          ENDDO
1518       ENDDO
1519!
1520!--    Count the total number of holes, required for informative message.
1521#if defined( __parallel )
1522       CALL MPI_ALLREDUCE( num_hole_l, num_hole, 1, MPI_INTEGER, MPI_SUM,      &
1523                           comm2d, ierr )
1524#else
1525       num_hole = num_hole_l
1526#endif   
1527       IF ( num_hole > 0  .AND.  .NOT. filled )  filled = .TRUE.
1528
1529    ENDDO
1530!
1531!-- Create an informative message if any holes were filled.
1532    IF ( filled )  THEN
1533       WRITE( message_string, * ) 'Topography was filtered, i.e. holes ' //    &
1534                                  'resolved by only one grid point '     //    &
1535                                  'were filled during initialization.'
1536       CALL message( 'init_grid', 'PA0430', 0, 0, 0, 6, 0 )
1537    ENDIF
1538
1539    DEALLOCATE( topo_tmp )
1540!
1541!-- Finally, exchange topo_3d array again and if necessary set Neumann boundary
1542!-- condition in case of non-cyclic lateral boundaries.
1543    CALL exchange_horiz_int( topo_3d, nys, nyn, nxl, nxr, nzt, nbgp )
1544
1545    IF ( .NOT. bc_ns_cyc )  THEN
1546       IF ( nys == 0  )  topo_3d(:,-1,:)   = topo_3d(:,0,:)
1547       IF ( nyn == ny )  topo_3d(:,ny+1,:) = topo_3d(:,ny,:)
1548    ENDIF
1549
1550    IF ( .NOT. bc_lr_cyc )  THEN
1551       IF ( nxl == 0  )  topo_3d(:,:,-1)   = topo_3d(:,:,0)
1552       IF ( nxr == nx )  topo_3d(:,:,nx+1) = topo_3d(:,:,nx)         
1553    ENDIF
1554!
1555!-- Exchange building ID and type. Note, building_type is an one-byte variable.
1556    IF ( building_id_f%from_file )                                             &
1557       CALL exchange_horiz_2d_int( building_id_f%var, nys, nyn, nxl, nxr, nbgp )
1558    IF ( building_type_f%from_file )                                           &
1559       CALL exchange_horiz_2d_byte( building_type_f%var, nys, nyn, nxl, nxr, nbgp )
1560
1561 END SUBROUTINE filter_topography
1562
1563
1564! Description:
1565! -----------------------------------------------------------------------------!
1566!> Reads topography information from file or sets generic topography. Moreover,
1567!> all topography-relevant topography arrays are initialized, and grid flags
1568!> are set. 
1569!------------------------------------------------------------------------------!
1570 SUBROUTINE init_topo( topo )
1571
1572    USE arrays_3d,                                                             &
1573        ONLY:  zw
1574       
1575    USE control_parameters,                                                    &
1576        ONLY:  bc_lr_cyc, bc_ns_cyc, building_height, building_length_x,       &
1577               building_length_y, building_wall_left, building_wall_south,     &
1578               canyon_height, canyon_wall_left, canyon_wall_south,             &
1579               canyon_width_x, canyon_width_y, dp_level_ind_b, dz,             &
1580               message_string, topography, topography_grid_convention,         &
1581               tunnel_height, tunnel_length, tunnel_width_x, tunnel_width_y,   &
1582               tunnel_wall_depth
1583         
1584    USE grid_variables,                                                        &
1585        ONLY:  dx, dy
1586       
1587    USE indices,                                                               &
1588        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nz,   &
1589               nzb, nzt
1590   
1591    USE kinds
1592   
1593    USE netcdf_data_input_mod,                                                 &
1594        ONLY:  buildings_f, terrain_height_f 
1595
1596    USE pegrid
1597
1598    IMPLICIT NONE
1599
1600    INTEGER(iwp) ::  bh                !< temporary vertical index of building height
1601    INTEGER(iwp) ::  ngp_bx            !< grid point number of building size along x
1602    INTEGER(iwp) ::  ngp_by            !< grid point number of building size along y
1603    INTEGER(iwp) ::  index_left_bwall  !< index for left building wall
1604    INTEGER(iwp) ::  index_right_bwall !< index for right building wall
1605    INTEGER(iwp) ::  index_north_bwall !< index for north building wall
1606    INTEGER(iwp) ::  index_south_bwall !< index for south building wall
1607    INTEGER(iwp) ::  ch                !< temporary vertical index for canyon height
1608    INTEGER(iwp) ::  ngp_cx            !< grid point number of canyon size along x
1609    INTEGER(iwp) ::  ngp_cy            !< grid point number of canyon size along y
1610    INTEGER(iwp) ::  index_left_cwall  !< index for left canyon wall
1611    INTEGER(iwp) ::  index_right_cwall !< index for right canyon wall
1612    INTEGER(iwp) ::  index_north_cwall !< index for north canyon wall
1613    INTEGER(iwp) ::  index_south_cwall !< index for south canyon wall
1614    INTEGER(iwp) ::  i                 !< index variable along x
1615    INTEGER(iwp) ::  j                 !< index variable along y
1616    INTEGER(iwp) ::  k                 !< index variable along z
1617    INTEGER(iwp) ::  hv_in             !< heavyside function to model inner tunnel surface
1618    INTEGER(iwp) ::  hv_out            !< heavyside function to model outer tunnel surface
1619    INTEGER(iwp) ::  txe_out           !< end position of outer tunnel wall in x
1620    INTEGER(iwp) ::  txs_out           !< start position of outer tunnel wall in x
1621    INTEGER(iwp) ::  tye_out           !< end position of outer tunnel wall in y
1622    INTEGER(iwp) ::  tys_out           !< start position of outer tunnel wall in y
1623    INTEGER(iwp) ::  txe_in            !< end position of inner tunnel wall in x
1624    INTEGER(iwp) ::  txs_in            !< start position of inner tunnel wall in x
1625    INTEGER(iwp) ::  tye_in            !< end position of inner tunnel wall in y
1626    INTEGER(iwp) ::  tys_in            !< start position of inner tunnel wall in y
1627    INTEGER(iwp) ::  td                !< tunnel wall depth
1628    INTEGER(iwp) ::  th                !< height of outer tunnel wall
1629
1630    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nzb_local         !< index for topography top at cell-center
1631    INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  topo !< input array for 3D topography and dummy array for setting "outer"-flags
1632!
1633!-- Check for correct setting of the namelist parameter topography. If
1634!-- topography information is read from file but topography = 'flat',
1635!-- initialization does not work properly.
1636    IF ( ( buildings_f%from_file  .OR.  terrain_height_f%from_file )  .AND.    &
1637           TRIM( topography ) /= 'read_from_file' )  THEN
1638       message_string =  'If topography information is provided (via ' //      &
1639                         'Netcdf or ASCII input) topography = '        //      &
1640                         '"read_from_file" is required.'
1641       CALL message( 'init_grid', 'PA0437', 1, 2, 0, 6, 0 )     
1642    ENDIF
1643!
1644!-- Set outer and inner index arrays for non-flat topography.
1645!-- Here consistency checks concerning domain size and periodicity are
1646!-- necessary.
1647!-- Within this SELECT CASE structure only nzb_local is initialized
1648!-- individually depending on the chosen topography type, all other index
1649!-- arrays are initialized further below.
1650    SELECT CASE ( TRIM( topography ) )
1651
1652       CASE ( 'flat' )
1653!   
1654!--       Initialilize 3D topography array, used later for initializing flags
1655          topo(nzb+1:nzt+1,:,:) = IBSET( topo(nzb+1:nzt+1,:,:), 0 )
1656         
1657       CASE ( 'closed_channel' )
1658!   
1659!--       Initialilize 3D topography array, used later for initializing flags
1660          topo(nzb+1:nzt,:,:) = IBSET( topo(nzb+1:nzt,:,:), 0 ) 
1661
1662       CASE ( 'single_building' )
1663!
1664!--       Single rectangular building, by default centered in the middle of the
1665!--       total domain
1666          ngp_bx = NINT( building_length_x / dx )
1667          ngp_by = NINT( building_length_y / dy )
1668          bh  = MINLOC( ABS( zw - building_height ), 1 ) - 1
1669          IF ( ABS( zw(bh)   - building_height ) == &
1670               ABS( zw(bh+1) - building_height )    )  bh = bh + 1
1671          IF ( building_wall_left == 9999999.9_wp )  THEN
1672             building_wall_left = ( nx + 1 - ngp_bx ) / 2 * dx
1673          ENDIF
1674          index_left_bwall = NINT( building_wall_left / dx )
1675          index_right_bwall = index_left_bwall + ngp_bx
1676
1677          IF ( building_wall_south == 9999999.9_wp )  THEN
1678              building_wall_south = ( ny + 1 - ngp_by ) / 2 * dy
1679          ENDIF
1680          index_south_bwall = NINT( building_wall_south / dy )
1681          index_north_bwall = index_south_bwall + ngp_by
1682
1683!
1684!--       Building size has to meet some requirements
1685          IF ( ( index_left_bwall < 1 ) .OR. ( index_right_bwall > nx-1 ) .OR. &
1686               ( index_right_bwall < index_left_bwall+3 ) .OR.                 &
1687               ( index_south_bwall < 1 ) .OR. ( index_north_bwall > ny-1 ) .OR.&
1688               ( index_north_bwall < index_south_bwall+3 ) )  THEN
1689             WRITE( message_string, * ) 'inconsistent building parameters:',   &
1690                                      '&index_left_bwall=', index_left_bwall,  &
1691                                      'index_right_bwall=', index_right_bwall, &
1692                                      'index_south_bwall=', index_south_bwall, &
1693                                      'index_north_bwall=', index_north_bwall, &
1694                                      'nx=', nx, 'ny=', ny
1695             CALL message( 'init_grid', 'PA0203', 1, 2, 0, 6, 0 )
1696          ENDIF
1697
1698          ALLOCATE( nzb_local(nysg:nyng,nxlg:nxrg) )
1699          nzb_local = 0
1700!
1701!--       Define the building.
1702          IF ( index_left_bwall <= nxr  .AND.  index_right_bwall >= nxl  .AND. &
1703               index_south_bwall <= nyn  .AND.  index_north_bwall >= nys )     & 
1704             nzb_local(MAX(nys,index_south_bwall):MIN(nyn,index_north_bwall),  &
1705                       MAX(nxl,index_left_bwall):MIN(nxr,index_right_bwall)) = bh
1706!
1707!--       Set bit array on basis of nzb_local
1708          DO  i = nxl, nxr
1709             DO  j = nys, nyn
1710                topo(nzb_local(j,i)+1:nzt+1,j,i) =                             &
1711                                 IBSET( topo(nzb_local(j,i)+1:nzt+1,j,i), 0 ) 
1712             ENDDO
1713          ENDDO
1714       
1715          DEALLOCATE( nzb_local )
1716
1717          CALL exchange_horiz_int( topo, nys, nyn, nxl, nxr, nzt, nbgp )
1718!
1719!--       Set boundary conditions also for flags. Can be interpreted as Neumann
1720!--       boundary conditions for topography.
1721          IF ( .NOT. bc_ns_cyc )  THEN
1722             IF ( nys == 0  )  THEN
1723                DO  i = 1, nbgp     
1724                   topo(:,nys-i,:)   = topo(:,nys,:)
1725                ENDDO
1726             ENDIF
1727             IF ( nyn == ny )  THEN
1728                DO  i = 1, nbgp 
1729                   topo(:,nyn+i,:) = topo(:,nyn,:)
1730                ENDDO
1731             ENDIF
1732          ENDIF
1733          IF ( .NOT. bc_lr_cyc )  THEN
1734             IF ( nxl == 0  )  THEN
1735                DO  i = 1, nbgp   
1736                   topo(:,:,nxl-i)   = topo(:,:,nxl)
1737                ENDDO
1738             ENDIF
1739             IF ( nxr == nx )  THEN
1740                DO  i = 1, nbgp   
1741                   topo(:,:,nxr+i) = topo(:,:,nxr)     
1742                ENDDO
1743             ENDIF     
1744          ENDIF
1745
1746       CASE ( 'single_street_canyon' )
1747!
1748!--       Single quasi-2D street canyon of infinite length in x or y direction.
1749!--       The canyon is centered in the other direction by default.
1750          IF ( canyon_width_x /= 9999999.9_wp )  THEN
1751!
1752!--          Street canyon in y direction
1753             ngp_cx = NINT( canyon_width_x / dx )
1754             IF ( canyon_wall_left == 9999999.9_wp )  THEN
1755                canyon_wall_left = ( nx + 1 - ngp_cx ) / 2 * dx
1756             ENDIF
1757             index_left_cwall= NINT( canyon_wall_left / dx )
1758             index_right_cwall= index_left_cwall+ ngp_cx
1759          ELSEIF ( canyon_width_y /= 9999999.9_wp )  THEN
1760!
1761!--          Street canyon in x direction
1762             ngp_cy = NINT( canyon_width_y / dy )
1763             IF ( canyon_wall_south == 9999999.9_wp )  THEN
1764                canyon_wall_south = ( ny + 1 - ngp_cy ) / 2 * dy
1765             ENDIF
1766             index_south_cwall = NINT( canyon_wall_south / dy )
1767             index_north_cwall = index_south_cwall + ngp_cy
1768     
1769          ELSE
1770             
1771             message_string = 'no street canyon width given'
1772             CALL message( 'init_grid', 'PA0204', 1, 2, 0, 6, 0 )
1773 
1774          ENDIF
1775
1776          ch  = MINLOC( ABS( zw - canyon_height ), 1 ) - 1
1777          IF ( ABS( zw(ch)   - canyon_height ) == &
1778               ABS( zw(ch+1) - canyon_height )    )  ch = ch + 1
1779          dp_level_ind_b = ch
1780!
1781!--       Street canyon size has to meet some requirements
1782          IF ( canyon_width_x /= 9999999.9_wp )  THEN
1783             IF ( ( index_left_cwall< 1 ) .OR. ( index_right_cwall> nx-1 ) .OR.&
1784                  ( ngp_cx < 3 ) .OR. ( ch < 3 ) )  THEN
1785                WRITE( message_string, * ) 'inconsistent canyon parameters:',  &
1786                                           '&index_left_cwall=', index_left_cwall, &
1787                                           ' index_right_cwall=', index_right_cwall, &
1788                                           ' ngp_cx=', ngp_cx,                 &
1789                                           ' ch=', ch, ' nx=', nx, ' ny=', ny
1790                CALL message( 'init_grid', 'PA0205', 1, 2, 0, 6, 0 ) 
1791             ENDIF
1792          ELSEIF ( canyon_width_y /= 9999999.9_wp )  THEN
1793             IF ( ( index_south_cwall < 1 ) .OR.                               &
1794                  ( index_north_cwall > ny-1 ) .OR. ( ngp_cy < 3 ) .OR.        &
1795                  ( ch < 3 ) )  THEN
1796                WRITE( message_string, * ) 'inconsistent canyon parameters:',  &
1797                                           '&index_south_cwall=', index_south_cwall, & 
1798                                           ' index_north_cwall=', index_north_cwall, &
1799                                           ' ngp_cy=', ngp_cy,                 &
1800                                           ' ch=', ch, ' nx=', nx, ' ny=', ny
1801                CALL message( 'init_grid', 'PA0206', 1, 2, 0, 6, 0 ) 
1802             ENDIF
1803          ENDIF
1804          IF ( canyon_width_x /= 9999999.9_wp .AND.                            &                 
1805               canyon_width_y /= 9999999.9_wp )  THEN
1806             message_string = 'inconsistent canyon parameters:' //             &   
1807                              '&street canyon can only be oriented' //         &
1808                              ' either in x- or in y-direction'
1809             CALL message( 'init_grid', 'PA0207', 1, 2, 0, 6, 0 )
1810          ENDIF
1811
1812          ALLOCATE( nzb_local(nysg:nyng,nxlg:nxrg) )
1813          nzb_local = ch
1814          IF ( canyon_width_x /= 9999999.9_wp )  THEN
1815             IF ( index_left_cwall<= nxr  .AND.  index_right_cwall>= nxl )     &
1816                nzb_local(:,MAX(nxl,index_left_cwall+1):MIN(nxr,index_right_cwall-1)) = 0
1817          ELSEIF ( canyon_width_y /= 9999999.9_wp )  THEN
1818             IF ( index_south_cwall <= nyn  .AND.  index_north_cwall >= nys )  &         
1819                nzb_local(MAX(nys,index_south_cwall+1):MIN(nyn,index_north_cwall-1),:) = 0
1820          ENDIF
1821!
1822!--       Set bit array on basis of nzb_local
1823          DO  i = nxl, nxr
1824             DO  j = nys, nyn
1825                topo(nzb_local(j,i)+1:nzt+1,j,i) =                             &
1826                                 IBSET( topo(nzb_local(j,i)+1:nzt+1,j,i), 0 ) 
1827             ENDDO
1828          ENDDO
1829          DEALLOCATE( nzb_local )
1830
1831          CALL exchange_horiz_int( topo, nys, nyn, nxl, nxr, nzt, nbgp )
1832!
1833!--       Set boundary conditions also for flags. Can be interpreted as Neumann
1834!--       boundary conditions for topography.
1835          IF ( .NOT. bc_ns_cyc )  THEN
1836             IF ( nys == 0  )  THEN
1837                DO  i = 1, nbgp     
1838                   topo(:,nys-i,:)   = topo(:,nys,:)
1839                ENDDO
1840             ENDIF
1841             IF ( nyn == ny )  THEN
1842                DO  i = 1, nbgp 
1843                   topo(:,nyn+i,:) = topo(:,nyn,:)
1844                ENDDO
1845             ENDIF
1846          ENDIF
1847          IF ( .NOT. bc_lr_cyc )  THEN
1848             IF ( nxl == 0  )  THEN
1849                DO  i = 1, nbgp   
1850                   topo(:,:,nxl-i)   = topo(:,:,nxl)
1851                ENDDO
1852             ENDIF
1853             IF ( nxr == nx )  THEN
1854                DO  i = 1, nbgp   
1855                   topo(:,:,nxr+i) = topo(:,:,nxr)     
1856                ENDDO
1857             ENDIF     
1858          ENDIF
1859
1860       CASE ( 'tunnel' )
1861
1862!
1863!--       Tunnel height
1864          IF ( tunnel_height == 9999999.9_wp )  THEN
1865             th = zw( INT( 0.2 * nz) )
1866          ELSE
1867             th = tunnel_height
1868          ENDIF
1869!
1870!--       Tunnel-wall depth
1871          IF ( tunnel_wall_depth == 9999999.9_wp )  THEN 
1872             td = MAX ( dx, dy, dz(1) )
1873          ELSE
1874             td = tunnel_wall_depth
1875          ENDIF
1876!
1877!--       Check for tunnel width
1878          IF ( tunnel_width_x == 9999999.9_wp  .AND.                           &
1879               tunnel_width_y == 9999999.9_wp  )  THEN
1880             message_string = 'No tunnel width is given. '
1881             CALL message( 'init_grid', 'PA0280', 1, 2, 0, 6, 0 )
1882          ENDIF
1883          IF ( tunnel_width_x /= 9999999.9_wp  .AND.                           &
1884               tunnel_width_y /= 9999999.9_wp  )  THEN
1885             message_string = 'Inconsistent tunnel parameters:' //             &   
1886                              'tunnel can only be oriented' //                 &
1887                              'either in x- or in y-direction.'
1888             CALL message( 'init_grid', 'PA0281', 1, 2, 0, 6, 0 )
1889          ENDIF
1890!
1891!--       Check for too small tunnel width in x- and y-direction
1892          IF ( tunnel_width_x /= 9999999.9_wp  .AND.                           &   
1893               tunnel_width_x - 2.0_wp * td <= 2.0_wp * dx )  THEN
1894             message_string = 'tunnel_width_x too small'
1895             CALL message( 'init_grid', 'PA0175', 1, 2, 0, 6, 0 )
1896          ENDIF
1897          IF ( tunnel_width_y /= 9999999.9_wp  .AND.                           &
1898               tunnel_width_y - 2.0_wp * td <= 2.0_wp * dy )  THEN
1899             message_string = 'tunnel_width_y too small'
1900             CALL message( 'init_grid', 'PA0455', 1, 2, 0, 6, 0 )
1901          ENDIF
1902!
1903!--       Check for too large tunnel width.
1904!--       Tunnel axis along y.
1905          IF ( tunnel_width_x /= 9999999.9_wp )  THEN
1906             IF ( tunnel_width_x > ( nx + 1 ) * dx )  THEN
1907                message_string = 'tunnel_width_x too large'
1908                CALL message( 'init_grid', 'PA0282', 1, 2, 0, 6, 0 )
1909             ENDIF
1910
1911             txs_out = INT( ( nx + 1 ) * 0.5_wp * dx - tunnel_width_x * 0.5_wp )
1912             txe_out = INT( ( nx + 1 ) * 0.5_wp * dx + tunnel_width_x * 0.5_wp )
1913             txs_in  = INT( ( nx + 1 ) * 0.5_wp * dx -                         &
1914                                      ( tunnel_width_x * 0.5_wp - td ) )
1915             txe_in  = INT( ( nx + 1 ) * 0.5_wp * dx +                         &
1916                                   ( tunnel_width_x * 0.5_wp - td ) )
1917
1918             tys_out = INT( ( ny + 1 ) * 0.5_wp * dy - tunnel_length * 0.5_wp )
1919             tye_out = INT( ( ny + 1 ) * 0.5_wp * dy + tunnel_length * 0.5_wp )
1920             tys_in  = tys_out
1921             tye_in  = tye_out
1922          ENDIF
1923!
1924!--       Tunnel axis along x.
1925          IF ( tunnel_width_y /= 9999999.9_wp )  THEN
1926             IF ( tunnel_width_y > ( ny + 1 ) * dy )  THEN
1927                message_string = 'tunnel_width_y too large'
1928                CALL message( 'init_grid', 'PA0456', 1, 2, 0, 6, 0 )
1929             ENDIF
1930
1931             txs_out = INT( ( nx + 1 ) * 0.5_wp * dx - tunnel_length * 0.5_wp )
1932             txe_out = INT( ( nx + 1 ) * 0.5_wp * dx + tunnel_length * 0.5_wp )
1933             txs_in  = txs_out
1934             txe_in  = txe_out
1935
1936             tys_out = INT( ( ny + 1 ) * 0.5_wp * dy - tunnel_width_y * 0.5_wp )
1937             tye_out = INT( ( ny + 1 ) * 0.5_wp * dy + tunnel_width_y * 0.5_wp )
1938             tys_in  = INT( ( ny + 1 ) * 0.5_wp * dy -                         &
1939                                        ( tunnel_width_y * 0.5_wp - td ) )
1940             tye_in  = INT( ( ny + 1 ) * 0.5_wp * dy +                         &
1941                                     ( tunnel_width_y * 0.5_wp - td ) )
1942          ENDIF
1943
1944          topo = 0
1945          DO  i = nxl, nxr
1946             DO  j = nys, nyn
1947!
1948!--             Use heaviside function to model outer tunnel surface
1949                hv_out = th * 0.5_wp *                                         &
1950                              ( ( SIGN( 1.0_wp, i * dx - txs_out ) + 1.0_wp )  &
1951                              - ( SIGN( 1.0_wp, i * dx - txe_out ) + 1.0_wp ) )
1952
1953                hv_out = hv_out * 0.5_wp *                                     &
1954                            ( ( SIGN( 1.0_wp, j * dy - tys_out ) + 1.0_wp )    &
1955                            - ( SIGN( 1.0_wp, j * dy - tye_out ) + 1.0_wp ) )
1956!   
1957!--             Use heaviside function to model inner tunnel surface
1958                hv_in  = ( th - td ) * 0.5_wp *                                &
1959                                ( ( SIGN( 1.0_wp, i * dx - txs_in ) + 1.0_wp ) &
1960                                - ( SIGN( 1.0_wp, i * dx - txe_in ) + 1.0_wp ) )
1961
1962                hv_in = hv_in * 0.5_wp *                                       &
1963                                ( ( SIGN( 1.0_wp, j * dy - tys_in ) + 1.0_wp ) &
1964                                - ( SIGN( 1.0_wp, j * dy - tye_in ) + 1.0_wp ) )
1965!
1966!--             Set flags at x-y-positions without any tunnel surface
1967                IF ( hv_out - hv_in == 0.0_wp )  THEN
1968                   topo(nzb+1:nzt+1,j,i) = IBSET( topo(nzb+1:nzt+1,j,i), 0 )
1969!
1970!--             Set flags at x-y-positions with tunnel surfaces
1971                ELSE
1972                   DO  k = nzb + 1, nzt + 1
1973!
1974!--                   Inner tunnel
1975                      IF ( hv_out - hv_in == th )  THEN
1976                         IF ( zw(k) <= hv_out )  THEN
1977                            topo(k,j,i) = IBCLR( topo(k,j,i), 0 )
1978                         ELSE
1979                            topo(k,j,i) = IBSET( topo(k,j,i), 0 )
1980                         ENDIF
1981                      ENDIF
1982!
1983!--                   Lateral tunnel walls
1984                      IF ( hv_out - hv_in == td )  THEN
1985                         IF ( zw(k) <= hv_in )  THEN
1986                            topo(k,j,i) = IBSET( topo(k,j,i), 0 )
1987                         ELSEIF ( zw(k) > hv_in  .AND.  zw(k) <= hv_out )  THEN
1988                            topo(k,j,i) = IBCLR( topo(k,j,i), 0 )
1989                         ELSEIF ( zw(k) > hv_out )  THEN
1990                            topo(k,j,i) = IBSET( topo(k,j,i), 0 )
1991                         ENDIF
1992                      ENDIF
1993                   ENDDO
1994                ENDIF
1995             ENDDO
1996          ENDDO
1997
1998          CALL exchange_horiz_int( topo, nys, nyn, nxl, nxr, nzt, nbgp )
1999!
2000!--       Set boundary conditions also for flags. Can be interpreted as Neumann
2001!--       boundary conditions for topography.
2002          IF ( .NOT. bc_ns_cyc )  THEN
2003             IF ( nys == 0  )  THEN
2004                DO  i = 1, nbgp     
2005                   topo(:,nys-i,:)   = topo(:,nys,:)
2006                ENDDO
2007             ENDIF
2008             IF ( nyn == ny )  THEN
2009                DO  i = 1, nbgp 
2010                   topo(:,nyn+i,:) = topo(:,nyn,:)
2011                ENDDO
2012             ENDIF
2013          ENDIF
2014          IF ( .NOT. bc_lr_cyc )  THEN
2015             IF ( nxl == 0  )  THEN
2016                DO  i = 1, nbgp   
2017                   topo(:,:,nxl-i)   = topo(:,:,nxl)
2018                ENDDO
2019             ENDIF
2020             IF ( nxr == nx )  THEN
2021                DO  i = 1, nbgp   
2022                   topo(:,:,nxr+i) = topo(:,:,nxr)     
2023                ENDDO
2024             ENDIF     
2025          ENDIF
2026
2027       CASE ( 'read_from_file' )
2028!
2029!--       Note, topography information have been already read. 
2030!--       If required, further process topography, i.e. reference buildings on
2031!--       top of orography and set temporary 3D topography array, which is
2032!--       used later to set grid flags. Calling of this rouinte is also
2033!--       required in case of ASCII input, even though no distinction between
2034!--       terrain- and building height is made in this case. 
2035          CALL process_topography( topo )
2036!
2037!--       Filter holes resolved by only one grid-point
2038          CALL filter_topography( topo )
2039!
2040!--       Exchange ghost-points, as well as add cyclic or Neumann boundary
2041!--       conditions.
2042          CALL exchange_horiz_int( topo, nys, nyn, nxl, nxr, nzt, nbgp )
2043!
2044!--       Set lateral boundary conditions for topography on all ghost layers
2045          IF ( .NOT. bc_ns_cyc )  THEN
2046             IF ( nys == 0  )  THEN
2047                DO  i = 1, nbgp         
2048                   topo(:,nys-i,:) = topo(:,nys,:)
2049                ENDDO
2050             ENDIF
2051             IF ( nyn == ny )  THEN
2052                DO  i = 1, nbgp         
2053                   topo(:,nyn+i,:) = topo(:,nyn,:)
2054                ENDDO
2055             ENDIF
2056          ENDIF
2057
2058          IF ( .NOT. bc_lr_cyc )  THEN
2059             IF ( nxl == 0  )  THEN
2060                DO  i = 1, nbgp 
2061                   topo(:,:,nxl-i) = topo(:,:,nxl)
2062                ENDDO
2063             ENDIF
2064             IF ( nxr == nx )  THEN
2065                DO  i = 1, nbgp 
2066                   topo(:,:,nxr+i) = topo(:,:,nxr)
2067                ENDDO
2068             ENDIF
2069          ENDIF
2070
2071
2072       CASE DEFAULT
2073!   
2074!--       The DEFAULT case is reached either if the parameter topography
2075!--       contains a wrong character string or if the user has defined a special
2076!--       case in the user interface. There, the subroutine user_init_grid
2077!--       checks which of these two conditions applies.
2078          CALL user_init_grid( topo )
2079          CALL filter_topography( topo )
2080
2081    END SELECT
2082!
2083!-- Consistency checks and index array initialization are only required for
2084!-- non-flat topography.
2085    IF ( TRIM( topography ) /= 'flat' )  THEN
2086!
2087!--    In case of non-flat topography, check whether the convention how to
2088!--    define the topography grid has been set correctly, or whether the default
2089!--    is applicable. If this is not possible, abort.
2090       IF ( TRIM( topography_grid_convention ) == ' ' )  THEN
2091          IF ( TRIM( topography ) /= 'closed_channel' .AND.                    &
2092               TRIM( topography ) /= 'single_building' .AND.                   &
2093               TRIM( topography ) /= 'single_street_canyon' .AND.              &
2094               TRIM( topography ) /= 'tunnel'  .AND.                           &
2095               TRIM( topography ) /= 'read_from_file')  THEN
2096!--          The default value is not applicable here, because it is only valid
2097!--          for the four standard cases 'single_building',
2098!--          'single_street_canyon', 'tunnel' and 'read_from_file'
2099!--          defined in init_grid.
2100             WRITE( message_string, * )                                        &
2101               'The value for "topography_grid_convention" ',                  &
2102               'is not set. Its default value is & only valid for ',           &
2103               '"topography" = ''single_building'', ''tunnel'' ',              &
2104               '''single_street_canyon'', ''closed_channel'' & or ',           &
2105               '''read_from_file''.',                                          &
2106               '& Choose ''cell_edge'' or ''cell_center''.'
2107             CALL message( 'init_grid', 'PA0239', 1, 2, 0, 6, 0 )
2108          ELSE
2109!--          The default value is applicable here.
2110!--          Set convention according to topography.
2111             IF ( TRIM( topography ) == 'single_building' .OR.                 &
2112                  TRIM( topography ) == 'single_street_canyon' )  THEN
2113                topography_grid_convention = 'cell_edge'
2114             ELSEIF ( TRIM( topography ) == 'read_from_file'  .OR.             &
2115                      TRIM( topography ) == 'tunnel')  THEN
2116                topography_grid_convention = 'cell_center'
2117             ENDIF
2118          ENDIF
2119       ELSEIF ( TRIM( topography_grid_convention ) /= 'cell_edge' .AND.        &
2120                TRIM( topography_grid_convention ) /= 'cell_center' )  THEN
2121          WRITE( message_string, * )                                           &
2122            'The value for "topography_grid_convention" is ',                  &
2123            'not recognized.& Choose ''cell_edge'' or ''cell_center''.'
2124          CALL message( 'init_grid', 'PA0240', 1, 2, 0, 6, 0 )
2125       ENDIF
2126
2127
2128       IF ( topography_grid_convention == 'cell_edge' )  THEN
2129!
2130!--       The array nzb_local as defined using the 'cell_edge' convention
2131!--       describes the actual total size of topography which is defined at the
2132!--       cell edges where u=0 on the topography walls in x-direction and v=0
2133!--       on the topography walls in y-direction. However, PALM uses individual
2134!--       arrays nzb_u|v|w|s_inner|outer that are based on nzb_s_inner.
2135!--       Therefore, the extent of topography in nzb_local is now reduced by
2136!--       1dx at the E topography walls and by 1dy at the N topography walls
2137!--       to form the basis for nzb_s_inner.
2138!--       Note, the reverse memory access (i-j instead of j-i) is absolutely
2139!--       required at this point.
2140          DO  j = nys+1, nyn+1
2141             DO  i = nxl-1, nxr
2142                DO  k = nzb, nzt+1
2143                   IF ( BTEST( topo(k,j,i), 0 )  .OR.                          &
2144                        BTEST( topo(k,j,i+1), 0 ) )                            &
2145                       topo(k,j,i) = IBSET( topo(k,j,i), 0 )
2146                ENDDO
2147             ENDDO
2148          ENDDO     
2149          CALL exchange_horiz_int( topo, nys, nyn, nxl, nxr, nzt, nbgp )
2150
2151          DO  i = nxl, nxr+1
2152             DO  j = nys-1, nyn
2153                DO  k = nzb, nzt+1
2154                   IF ( BTEST( topo(k,j,i), 0 )  .OR.                          &
2155                        BTEST( topo(k,j+1,i), 0 ) )                            &
2156                      topo(k,j,i) = IBSET( topo(k,j,i), 0 )
2157                ENDDO
2158             ENDDO
2159          ENDDO 
2160          CALL exchange_horiz_int( topo, nys, nyn, nxl, nxr, nzt, nbgp )
2161   
2162       ENDIF
2163    ENDIF
2164
2165
2166 END SUBROUTINE init_topo
2167
2168 SUBROUTINE set_topo_flags(topo)
2169
2170    USE control_parameters,                                                    &
2171        ONLY:  bc_lr_cyc, bc_ns_cyc, constant_flux_layer,                      &
2172               scalar_advec, topography, use_surface_fluxes, use_top_fluxes
2173
2174    USE indices,                                                               &
2175        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb,  &
2176               nzt, topo_top_ind, wall_flags_static_0, wall_flags_total_0
2177
2178    USE kinds
2179
2180    IMPLICIT NONE
2181
2182    INTEGER(iwp) ::  i             !< index variable along x
2183    INTEGER(iwp) ::  ibit          !< integer bit position of topgraphy masking array
2184    INTEGER(iwp) ::  j             !< index variable along y
2185    INTEGER(iwp) ::  k             !< index variable along z
2186
2187    INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  topo !< input array for 3D topography and dummy array for setting "outer"-flags
2188
2189    ALLOCATE( wall_flags_static_0(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2190    wall_flags_static_0 = 0
2191!
2192!-- Set-up topography flags. First, set flags only for s, u, v and w-grid.
2193!-- Further special flags will be set in following loops.
2194    DO  i = nxl, nxr
2195       DO  j = nys, nyn
2196          DO  k = nzb, nzt+1
2197!
2198!--          scalar grid
2199             IF ( BTEST( topo(k,j,i), 0 ) )                                    &
2200                wall_flags_static_0(k,j,i) = IBSET( wall_flags_static_0(k,j,i), 0 )
2201!
2202!--          u grid
2203             IF ( BTEST( topo(k,j,i),   0 )  .AND.                             &
2204                  BTEST( topo(k,j,i-1), 0 ) )                                  &
2205                wall_flags_static_0(k,j,i) = IBSET( wall_flags_static_0(k,j,i), 1 )
2206!
2207!--          v grid
2208             IF ( BTEST( topo(k,j,i),   0 )  .AND.                             &
2209                  BTEST( topo(k,j-1,i), 0 ) )                                  &
2210                 wall_flags_static_0(k,j,i) = IBSET( wall_flags_static_0(k,j,i), 2 )
2211
2212          ENDDO
2213
2214          DO k = nzb, nzt
2215!
2216!--          w grid
2217             IF ( BTEST( topo(k,j,i),   0 )  .AND.                             &
2218                  BTEST( topo(k+1,j,i), 0 ) )                                  &
2219                wall_flags_static_0(k,j,i) = IBSET( wall_flags_static_0(k,j,i), 3 )
2220          ENDDO
2221         
2222          IF ( topography /= 'closed_channel' ) THEN
2223             wall_flags_static_0(nzt+1,j,i) = IBSET( wall_flags_static_0(nzt+1,j,i), 3 )
2224          ENDIF
2225
2226       ENDDO
2227    ENDDO
2228
2229    CALL exchange_horiz_int( wall_flags_static_0, nys, nyn, nxl, nxr, nzt, nbgp )
2230
2231!
2232!-- Set outer array for scalars to mask near-surface grid points. Note, on
2233!-- basis of flag 24 futher flags will be derived which are used to control
2234!-- production of subgrid TKE production near walls.
2235   
2236    ALLOCATE( wall_flags_total_0(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2237    wall_flags_total_0 = 0
2238                                   
2239    DO i = nxl, nxr
2240       DO j = nys, nyn
2241          DO k = nzb, nzt+1
2242             wall_flags_total_0(k,j,i) = IOR( wall_flags_total_0(k,j,i), wall_flags_static_0(k,j,i) )
2243          ENDDO
2244       ENDDO
2245    ENDDO
2246   
2247    CALL exchange_horiz_int( wall_flags_total_0, nys, nyn, nxl, nxr, nzt, nbgp )
2248   
2249    DO i = nxl, nxr
2250       DO j = nys, nyn
2251          DO k = nzb, nzt+1
2252             IF ( BTEST( wall_flags_total_0(k,j-1,i), 0 )    .AND.                   &
2253                  BTEST( wall_flags_total_0(k,j+1,i), 0 )    .AND.                   &
2254                  BTEST( wall_flags_total_0(k,j,i-1), 0 )    .AND.                   &
2255                  BTEST( wall_flags_total_0(k,j,i+1), 0 )    .AND.                   &
2256                  BTEST( wall_flags_total_0(k,j-1,i-1), 0 )  .AND.                   &
2257                  BTEST( wall_flags_total_0(k,j+1,i-1), 0 )  .AND.                   &
2258                  BTEST( wall_flags_total_0(k,j-1,i+1), 0 )  .AND.                   &
2259                  BTEST( wall_flags_total_0(k,j+1,i+1), 0 ) )                        &
2260                wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 24 )
2261          ENDDO
2262       ENDDO
2263    ENDDO
2264!
2265!-- Set further special flags
2266    DO i = nxl, nxr
2267       DO j = nys, nyn
2268          DO k = nzb, nzt+1
2269!
2270!--          scalar grid, former nzb_diff_s_inner.
2271!--          Note, use this flag also to mask topography in diffusion_u and
2272!--          diffusion_v along the vertical direction. In case of
2273!--          use_surface_fluxes, fluxes are calculated via MOST, else, simple
2274!--          gradient approach is applied. Please note, in case of u- and v-
2275!--          diffuison, a small error is made at edges (on the east side for u,
2276!--          at the north side for v), since topography on scalar grid point
2277!--          is used instead of topography on u/v-grid. As number of topography grid
2278!--          points on uv-grid is different than s-grid, different number of
2279!--          surface elements would be required. In order to avoid this,
2280!--          treat edges (u(k,j,i+1)) simply by a gradient approach, i.e. these
2281!--          points are not masked within diffusion_u. Tests had shown that the
2282!--          effect on the flow is negligible.
2283             IF ( constant_flux_layer  .OR.  use_surface_fluxes )  THEN
2284                IF ( BTEST( wall_flags_total_0(k,j,i), 0 ) )                         &
2285                   wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 8 )
2286             ELSE
2287                wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 8 )
2288             ENDIF
2289
2290          ENDDO
2291!
2292!--       Special flag to control vertical diffusion at model top - former
2293!--       nzt_diff
2294          wall_flags_total_0(:,j,i) = IBSET( wall_flags_total_0(:,j,i), 9 )
2295          IF ( use_top_fluxes )                                                &
2296             wall_flags_total_0(nzt+1,j,i) = IBCLR( wall_flags_total_0(nzt+1,j,i), 9 )
2297
2298
2299          DO k = nzb+1, nzt
2300!
2301!--          Special flag on u grid, former nzb_u_inner + 1, required   
2302!--          for disturb_field and initialization. Do not disturb directly at
2303!--          topography, as well as initialize u with zero one grid point outside
2304!--          of topography.
2305             IF ( BTEST( wall_flags_total_0(k-1,j,i), 1 )  .AND.                     &
2306                  BTEST( wall_flags_total_0(k,j,i),   1 )  .AND.                     &
2307                  BTEST( wall_flags_total_0(k+1,j,i), 1 ) )                          &
2308                wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 20 )
2309!
2310!--          Special flag on v grid, former nzb_v_inner + 1, required   
2311!--          for disturb_field and initialization. Do not disturb directly at
2312!--          topography, as well as initialize v with zero one grid point outside
2313!--          of topography.
2314             IF ( BTEST( wall_flags_total_0(k-1,j,i), 2 )  .AND.                     &
2315                  BTEST( wall_flags_total_0(k,j,i),   2 )  .AND.                     &
2316                  BTEST( wall_flags_total_0(k+1,j,i), 2 ) )                          &
2317                wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 21 )
2318!
2319!--          Special flag on scalar grid, former nzb_s_inner+1. Used for
2320!--          lpm_sgs_tke
2321             IF ( BTEST( wall_flags_total_0(k,j,i),   0 )  .AND.                     &
2322                  BTEST( wall_flags_total_0(k-1,j,i), 0 )  .AND.                     &
2323                  BTEST( wall_flags_total_0(k+1,j,i), 0 ) )                          &
2324                wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 25 )
2325!
2326!--          Special flag on scalar grid, nzb_diff_s_outer - 1, required in
2327!--          in production_e
2328             IF ( constant_flux_layer  .OR.  use_surface_fluxes )  THEN
2329                IF ( BTEST( wall_flags_total_0(k,j,i),   24 )  .AND.                 &
2330                     BTEST( wall_flags_total_0(k-1,j,i), 24 )  .AND.                 &
2331                     BTEST( wall_flags_total_0(k+1,j,i), 0 ) )                       &
2332                   wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 29 )
2333             ELSE
2334                IF ( BTEST( wall_flags_total_0(k,j,i), 0 ) )                         &
2335                   wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 29 )
2336             ENDIF
2337!
2338!--          Special flag on scalar grid, nzb_diff_s_outer - 1, required in
2339!--          in production_e
2340             IF ( constant_flux_layer  .OR.  use_surface_fluxes )  THEN
2341                IF ( BTEST( wall_flags_total_0(k,j,i),   0 )  .AND.                  &
2342                     BTEST( wall_flags_total_0(k-1,j,i), 0 )  .AND.                  &
2343                     BTEST( wall_flags_total_0(k+1,j,i), 0 ) )                       &
2344                   wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 30 )
2345             ELSE
2346                IF ( BTEST( wall_flags_total_0(k,j,i), 0 ) )                         &
2347                   wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 30 )
2348             ENDIF
2349          ENDDO
2350!
2351!--       Flags indicating downward facing walls
2352          DO k = nzb+1, nzt+1
2353!
2354!--          Scalar grid
2355             IF ( BTEST( wall_flags_total_0(k-1,j,i), 0 )  .AND.                     &
2356            .NOT. BTEST( wall_flags_total_0(k,j,i), 0   ) )                          & 
2357                 wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 13 ) 
2358!
2359!--          Downward facing wall on u grid
2360             IF ( BTEST( wall_flags_total_0(k-1,j,i), 1 )  .AND.                     &
2361            .NOT. BTEST( wall_flags_total_0(k,j,i), 1   ) )                          & 
2362                wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 15 )
2363!
2364!--          Downward facing wall on v grid
2365             IF ( BTEST( wall_flags_total_0(k-1,j,i), 2 )  .AND.                     &
2366            .NOT. BTEST( wall_flags_total_0(k,j,i), 2   ) )                          & 
2367                wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 17 )
2368!
2369!--          Downward facing wall on w grid
2370             IF ( BTEST( wall_flags_total_0(k-1,j,i), 3 )  .AND.                     &
2371            .NOT. BTEST( wall_flags_total_0(k,j,i), 3 ) )                            & 
2372                wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 19 )
2373          ENDDO
2374!
2375!--       Flags indicating upward facing walls
2376          DO k = nzb, nzt
2377!
2378!--          Upward facing wall on scalar grid
2379             IF ( .NOT. BTEST( wall_flags_total_0(k,j,i),   0 )  .AND.               &
2380                        BTEST( wall_flags_total_0(k+1,j,i), 0 ) )                    & 
2381                wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 12 )
2382!
2383!--          Upward facing wall on u grid
2384             IF ( .NOT. BTEST( wall_flags_total_0(k,j,i),   1 )  .AND.               &
2385                        BTEST( wall_flags_total_0(k+1,j,i), 1 ) )                    & 
2386                wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 14 )
2387
2388!   
2389!--          Upward facing wall on v grid
2390             IF ( .NOT. BTEST( wall_flags_total_0(k,j,i),   2 )  .AND.               &
2391                        BTEST( wall_flags_total_0(k+1,j,i), 2 ) )                    & 
2392                wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 16 )
2393   
2394!
2395!--          Upward facing wall on w grid
2396             IF ( .NOT. BTEST( wall_flags_total_0(k,j,i),   3 )  .AND.               &
2397                        BTEST( wall_flags_total_0(k+1,j,i), 3 ) )                    & 
2398                wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 18 )
2399!
2400!--          Special flag on scalar grid, former nzb_s_inner
2401             IF ( BTEST( wall_flags_total_0(k,j,i), 0 )  .OR.                        &
2402                  BTEST( wall_flags_total_0(k,j,i), 12 ) .OR.                        &
2403                  BTEST( wall_flags_total_0(k,j,i), 13 ) )                           &
2404                   wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 22 )
2405!
2406!--          Special flag on scalar grid, nzb_diff_s_inner - 1, required for
2407!--          flow_statistics
2408             IF ( constant_flux_layer  .OR.  use_surface_fluxes )  THEN
2409                IF ( BTEST( wall_flags_total_0(k,j,i),   0 )  .AND.                  &
2410                     BTEST( wall_flags_total_0(k+1,j,i), 0 ) )                       &
2411                  wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 23 )
2412             ELSE
2413                IF ( BTEST( wall_flags_total_0(k,j,i), 22 ) )                        &
2414                   wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 23 )
2415             ENDIF
2416   
2417
2418          ENDDO
2419          wall_flags_total_0(nzt+1,j,i) = IBSET( wall_flags_total_0(nzt+1,j,i), 22 )
2420          wall_flags_total_0(nzt+1,j,i) = IBSET( wall_flags_total_0(nzt+1,j,i), 23 )
2421!
2422!--       Set flags indicating that topography is close by in horizontal
2423!--       direction, i.e. flags that infold the topography. These will be used
2424!--       to set advection flags for passive scalars, where due to large
2425!--       gradients near buildings stationary numerical oscillations can produce
2426!--       unrealistically high concentrations. This is only necessary if
2427!--       WS-scheme is applied for scalar advection. Note, these flags will be
2428!--       only used for passive scalars such as chemical species or aerosols.
2429          IF ( scalar_advec == 'ws-scheme' )  THEN
2430             DO k = nzb, nzt
2431                IF ( BTEST( wall_flags_total_0(k,j,i), 0 )  .AND. (                  &
2432                     ANY( .NOT. BTEST( wall_flags_total_0(k,j-3:j+3,i-1), 0 ) )  .OR.&
2433                     ANY( .NOT. BTEST( wall_flags_total_0(k,j-3:j+3,i-2), 0 ) )  .OR.&
2434                     ANY( .NOT. BTEST( wall_flags_total_0(k,j-3:j+3,i-3), 0 ) )  .OR.&
2435                     ANY( .NOT. BTEST( wall_flags_total_0(k,j-3:j+3,i+1), 0 ) )  .OR.&
2436                     ANY( .NOT. BTEST( wall_flags_total_0(k,j-3:j+3,i+2), 0 ) )  .OR.&
2437                     ANY( .NOT. BTEST( wall_flags_total_0(k,j-3:j+3,i+3), 0 ) )  .OR.&
2438                     ANY( .NOT. BTEST( wall_flags_total_0(k,j-1,i-3:i+3), 0 ) )  .OR.&
2439                     ANY( .NOT. BTEST( wall_flags_total_0(k,j-2,i-3:i+3), 0 ) )  .OR.&
2440                     ANY( .NOT. BTEST( wall_flags_total_0(k,j-3,i-3:i+3), 0 ) )  .OR.&
2441                     ANY( .NOT. BTEST( wall_flags_total_0(k,j+1,i-3:i+3), 0 ) )  .OR.&
2442                     ANY( .NOT. BTEST( wall_flags_total_0(k,j+2,i-3:i+3), 0 ) )  .OR.&
2443                     ANY( .NOT. BTEST( wall_flags_total_0(k,j+3,i-3:i+3), 0 ) )      &
2444                                                            ) )                      &
2445                   wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 31 )
2446                     
2447             ENDDO
2448          ENDIF
2449       ENDDO
2450    ENDDO
2451!
2452!-- Finally, set identification flags indicating natural terrain or buildings.
2453!-- Natural terrain grid points. Information on the type of the surface is
2454!-- stored in bit 1 of 3D Integer array topo. However, this bit is only set
2455!-- when topography is read from file. In order to run the land-surface model
2456!-- also without topography information, set bit 1 explicitely in this case.
2457!--
2458!-- Natural terrain grid points
2459!-- If no topography is initialized, the land-surface is at k = nzb.
2460    IF ( TRIM( topography ) /= 'read_from_file' )  THEN
2461       wall_flags_static_0(nzb,:,:) = IBSET( wall_flags_static_0(nzb,:,:), 5 )
2462    ELSE
2463       DO i = nxl, nxr
2464          DO j = nys, nyn
2465             DO k = nzb, nzt+1
2466!         
2467!--             Natural terrain grid point
2468                IF ( BTEST( topo(k,j,i), 1 ) )                                 &
2469                   wall_flags_static_0(k,j,i) = IBSET( wall_flags_static_0(k,j,i), 5 )
2470             ENDDO
2471          ENDDO
2472       ENDDO
2473    ENDIF
2474!
2475!-- Building grid points.
2476    DO i = nxl, nxr
2477       DO j = nys, nyn
2478          DO k = nzb, nzt+1
2479             IF ( BTEST( topo(k,j,i), 2 ) )                                    &
2480                wall_flags_static_0(k,j,i) = IBSET( wall_flags_static_0(k,j,i), 6 )
2481          ENDDO
2482       ENDDO
2483    ENDDO
2484!
2485!-- Set flag 4, indicating new topography grid points due to filtering.
2486    DO i = nxl, nxr
2487       DO j = nys, nyn
2488          DO k = nzb, nzt+1
2489             IF ( BTEST( topo(k,j,i), 4 ) )                                    &
2490                wall_flags_static_0(k,j,i) = IBSET( wall_flags_static_0(k,j,i), 4 )
2491          ENDDO
2492       ENDDO
2493    ENDDO
2494   
2495    CALL exchange_horiz_int( wall_flags_static_0, nys, nyn, nxl, nxr, nzt, nbgp )
2496   
2497    DO i = nxl, nxr
2498       DO j = nys, nyn
2499          DO k = nzb, nzt+1
2500             wall_flags_total_0(k,j,i) = IOR( wall_flags_total_0(k,j,i), wall_flags_static_0(k,j,i) )
2501          ENDDO
2502       ENDDO
2503    ENDDO
2504!
2505!-- Exchange ghost points for wall flags
2506    CALL exchange_horiz_int( wall_flags_total_0, nys, nyn, nxl, nxr, nzt, nbgp )
2507!
2508!-- Set boundary conditions also for flags. Can be interpreted as Neumann
2509!-- boundary conditions for topography.
2510    IF ( .NOT. bc_ns_cyc )  THEN
2511       IF ( nys == 0  )  THEN
2512          DO  i = 1, nbgp     
2513             wall_flags_total_0(:,nys-i,:)   = wall_flags_total_0(:,nys,:)
2514          ENDDO
2515       ENDIF
2516       IF ( nyn == ny )  THEN
2517          DO  i = 1, nbgp 
2518             wall_flags_total_0(:,nyn+i,:) = wall_flags_total_0(:,nyn,:)
2519          ENDDO
2520       ENDIF
2521    ENDIF
2522    IF ( .NOT. bc_lr_cyc )  THEN
2523       IF ( nxl == 0  )  THEN
2524          DO  i = 1, nbgp   
2525             wall_flags_total_0(:,:,nxl-i)   = wall_flags_total_0(:,:,nxl)
2526          ENDDO
2527       ENDIF
2528       IF ( nxr == nx )  THEN
2529          DO  i = 1, nbgp   
2530             wall_flags_total_0(:,:,nxr+i) = wall_flags_total_0(:,:,nxr)     
2531          ENDDO
2532       ENDIF     
2533    ENDIF
2534!
2535!-- Pre-calculate topography top indices (former get_topography_top_index
2536!-- function)
2537    ALLOCATE( topo_top_ind(nysg:nyng,nxlg:nxrg,0:4) )
2538!
2539!-- Uppermost topography index on scalar grid
2540    ibit = 12
2541    topo_top_ind(:,:,0) = MAXLOC(                                              &
2542                                  MERGE( 1, 0,                                 &
2543                                    BTEST( wall_flags_total_0(:,:,:), ibit )   &
2544                                       ), DIM = 1                              &
2545                                ) - 1 
2546!
2547!-- Uppermost topography index on u grid
2548    ibit = 14
2549    topo_top_ind(:,:,1) = MAXLOC(                                              &
2550                                  MERGE( 1, 0,                                 &
2551                                    BTEST( wall_flags_total_0(:,:,:), ibit )   &
2552                                       ), DIM = 1                              &
2553                                ) - 1 
2554!
2555!-- Uppermost topography index on v grid
2556    ibit = 16
2557    topo_top_ind(:,:,2) = MAXLOC(                                              &
2558                                  MERGE( 1, 0,                                 &
2559                                    BTEST( wall_flags_total_0(:,:,:), ibit )   &
2560                                       ), DIM = 1                              &
2561                                ) - 1 
2562!
2563!-- Uppermost topography index on w grid
2564    ibit = 18
2565    topo_top_ind(:,:,3) = MAXLOC(                                              &
2566                                  MERGE( 1, 0,                                 &
2567                                    BTEST( wall_flags_total_0(:,:,:), ibit )   &
2568                                       ), DIM = 1                              &
2569                                ) - 1 
2570!
2571!-- Uppermost topography index on scalar outer grid
2572    ibit = 24
2573    topo_top_ind(:,:,4) = MAXLOC(                                              &
2574                                  MERGE( 1, 0,                                 &
2575                                    BTEST( wall_flags_total_0(:,:,:), ibit )   &
2576                                       ), DIM = 1                              &
2577                                ) - 1
2578
2579 END SUBROUTINE set_topo_flags
2580
2581
2582
2583
Note: See TracBrowser for help on using the repository browser.