source: palm/trunk/SOURCE/poismg_noopt_mod.f90 @ 3182

Last change on this file since 3182 was 3182, checked in by suehring, 3 years ago

New Inifor features: grid stretching, improved command-interface, support start dates in different formats in both YYYYMMDD and YYYYMMDDHH, Ability to manually control input file prefixes (--radiation-prefix, --soil-preifx, --flow-prefix, --soilmoisture-prefix) for compatiblity with DWD forcast naming scheme; GNU-style short and long option; Prepared output of large-scale forcing profiles (no computation yet); Added preprocessor flag netcdf4 to switch output format between netCDF 3 and 4; Updated netCDF variable names and attributes to comply with PIDS v1.9; Inifor bugfixes: Improved compatibility with older Intel Intel compilers by avoiding implicit array allocation; Added origin_lon/_lat values and correct reference time in dynamic driver global attributes; corresponding PALM changes: adjustments to revised Inifor; variables names in dynamic driver adjusted; enable geostrophic forcing also in offline nested mode; variable names in LES-LES and COSMO offline nesting changed; lateral boundary flags for nesting, in- and outflow conditions renamed

  • Property svn:keywords set to Id
File size: 80.8 KB
Line 
1!> @file poismg_noopt_mod.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-2018 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22! Rename variables in mesoscale-offline nesting mode
23!
24! Former revisions:
25! -----------------
26! $Id: poismg_noopt_mod.f90 3182 2018-07-27 13:36:03Z suehring $
27! Set lateral boundary conditions for divergence
28!
29! 2793 2018-02-07 10:54:33Z suehring
30! Corrected "Former revisions" section
31!
32! 2696 2017-12-14 17:12:51Z kanani
33! - Change in file header (GPL part)
34! - Rename file into poismg_noopt_mod
35! - Modularize poismg_noopt - initialization of masking flags moved from
36!   init_grid
37! - Masking of downward-facing walls walls
38! - Large-scale forcing implemented (MS)
39!
40! 2591 2017-10-30 08:07:25Z maronga
41!
42! 2586 2017-10-27 13:03:50Z kanani
43!
44! 2232 2017-05-30 17:47:52Z suehring
45! Bugfixes OpenMP
46!
47! 2037 2016-10-26 11:15:40Z knoop
48! Anelastic approximation implemented
49!
50! 2021 2016-10-07 14:08:57Z suehring
51! Bugfix: restore nest_bound_(l/r/s/n) in case of mg_switch_to_pe0
52!
53! 2000 2016-08-20 18:09:15Z knoop
54! Forced header and separation lines into 80 columns
55!
56! 1934 2016-06-13 09:46:57Z hellstea
57! Rename subroutines and cpu-measure log points to indicate _noopt version
58!
59! 1762 2016-02-25 12:31:13Z hellstea
60! Introduction of nested domain feature
61!
62! 1682 2015-10-07 23:56:08Z knoop
63! Code annotations made doxygen readable
64!
65! 1353 2014-04-08 15:21:23Z heinze
66! REAL constants provided with KIND-attribute
67!
68! 1322 2014-03-20 16:38:49Z raasch
69! REAL constants defined as wp-kind
70!
71! 1320 2014-03-20 08:40:49Z raasch
72! ONLY-attribute added to USE-statements,
73! kind-parameters added to all INTEGER and REAL declaration statements,
74! kinds are defined in new module kinds,
75! old module precision_kind is removed,
76! revision history before 2012 removed,
77! comment fields (!:) to be used for variable explanations added to
78! all variable declaration statements
79!
80! 1318 2014-03-17 13:35:16Z raasch
81! module interfaces removed
82!
83! 1159 2013-05-21 11:58:22Z fricke
84! bc_lr/ns_dirneu/neudir removed
85!
86! 1092 2013-02-02 11:24:22Z raasch
87! unused variables removed
88!
89! 1056 2012-11-16 15:28:04Z raasch
90! Bugfix: all ghost points have to be used for allocating p3
91! arrays p2, f2, and f2_l changed from allocatable to automatic
92!
93! 1036 2012-10-22 13:43:42Z raasch
94! code put under GPL (PALM 3.9)
95!
96! 996 2012-09-07 10:41:47Z raasch
97! little reformatting
98!
99! 978 2012-08-09 08:28:32Z fricke
100! bc_lr/ns_dirneu/neudir added
101!
102! 880 2012-04-13 06:28:59Z raasch
103! Bugfix: preprocessor statements for parallel execution added
104!
105! 778 2011-11-07 14:18:25Z fricke
106! Allocation of p3 changes when multigrid is used and the collected field on PE0
107! has more grid points than the subdomain of an PE.
108!
109! 707 2011-03-29 11:39:40Z raasch
110! p_loc is used instead of p in the main routine (poismg).
111! On coarse grid levels, gathered data are identically processed on all PEs
112! (before, on PE0 only), so that the subsequent scattering of data is not
113! neccessary any more.
114! bc_lr/ns replaced by bc_lr/ns_cyc/dirrad/raddir
115! Bugfix: bottom (nzb) and top (nzt+1) boundary conditions set in routines
116! resid and restrict. They were missed before which may have led to
117! unpredictable results.
118!
119! 667 2010-12-23 12:06:00Z suehring/gryschka
120! Calls of exchange_horiz are modified.
121!
122! 622 2010-12-10 08:08:13Z raasch
123! optional barriers included in order to speed up collective operations
124!
125! 257 2009-03-11 15:17:42Z heinze
126! Output of messages replaced by message handling routine.
127!
128! 181 2008-07-30 07:07:47Z raasch
129! Bugfix: grid_level+1 has to be used in restrict for flags-array
130!
131! 114 2007-10-10 00:03:15Z raasch
132! Boundary conditions at walls are implicitly set using flag arrays. Only
133! Neumann BC is allowed. Upper walls are still not realized.
134! Bottom and top BCs for array f_mg in restrict removed because boundary
135! values are not needed (right hand side of SOR iteration).
136!
137! 75 2007-03-22 09:54:05Z raasch
138! 2nd+3rd argument removed from exchange horiz
139!
140! RCS Log replace by Id keyword, revision history cleaned up
141!
142! Revision 1.6  2005/03/26 20:55:54  raasch
143! Implementation of non-cyclic (Neumann) horizontal boundary conditions,
144! routine prolong simplified (one call of exchange_horiz spared)
145!
146! Revision 1.1  2001/07/20 13:10:51  raasch
147! Initial revision
148!
149!
150! Description:
151! ------------
152!> Solves the Poisson equation for the perturbation pressure with a multigrid
153!> V- or W-Cycle scheme.
154!>
155!> This multigrid method was originally developed for PALM by Joerg Uhlenbrock,
156!> September 2000 - July 2001.
157!>
158!> @attention Loop unrolling and cache optimization in SOR-Red/Black method
159!>            still does not give the expected speedup!
160!>
161!> @todo Further work required.
162!> @todo Formatting adjustments required (indention after modularization)
163!------------------------------------------------------------------------------!
164 MODULE poismg_noopt_mod
165 
166    USE control_parameters,                                                    &
167        ONLY:  bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r,                 &
168               bc_dirichlet_s, bc_radiation_l, bc_radiation_n, bc_radiation_r, &
169               bc_radiation_s, child_domain, grid_level, nesting_offline
170
171    USE cpulog,                                                                &
172        ONLY:  cpu_log, log_point_s
173
174    USE kinds
175
176    USE pegrid
177
178    PRIVATE
179
180    INTERFACE poismg_noopt
181       MODULE PROCEDURE poismg_noopt
182    END INTERFACE poismg_noopt
183
184    INTERFACE poismg_noopt_init
185       MODULE PROCEDURE poismg_noopt_init
186    END INTERFACE poismg_noopt_init
187
188    PUBLIC poismg_noopt, poismg_noopt_init
189
190 CONTAINS
191
192    SUBROUTINE poismg_noopt( r )
193   
194
195       USE arrays_3d,                                                          &
196           ONLY:  d, p_loc
197
198       USE control_parameters,                                                 &
199           ONLY:  bc_lr_cyc, bc_ns_cyc, gathered_size, grid_level_count,       &
200                  ibc_p_t, maximum_grid_level, message_string, mgcycles,       &
201                  mg_cycles, mg_switch_to_pe0_level, residual_limit,           &
202                  subdomain_size
203
204       USE control_parameters,                                                 &
205           ONLY:  bc_lr_cyc, bc_ns_cyc, gathered_size, grid_level,             &
206                  grid_level_count, ibc_p_t,                                   &
207                  maximum_grid_level, message_string, mgcycles, mg_cycles,     &
208                  mg_switch_to_pe0_level, residual_limit, subdomain_size
209
210       USE cpulog,                                                             &
211           ONLY:  cpu_log, log_point_s
212
213       USE indices,                                                            &
214           ONLY:  nxl, nxlg, nxl_mg, nxr, nxrg, nxr_mg, nys, nysg, nys_mg, nyn,&
215                  nyng, nyn_mg, nzb, nzt, nzt_mg
216
217       USE kinds
218
219       USE pegrid
220
221       IMPLICIT NONE
222
223       REAL(wp) ::  maxerror          !<
224       REAL(wp) ::  maximum_mgcycles  !<
225       REAL(wp) ::  residual_norm     !<
226
227       REAL(wp), DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) ::  r  !<
228
229       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  p3  !<
230
231
232       CALL cpu_log( log_point_s(29), 'poismg_noopt', 'start' )
233!
234!--    Initialize arrays and variables used in this subroutine
235
236!--    If the number of grid points of the gathered grid, which is collected
237!--    on PE0, is larger than the number of grid points of an PE, than array
238!--    p3 will be enlarged.
239       IF ( gathered_size > subdomain_size )  THEN
240          ALLOCATE( p3(nzb:nzt_mg(mg_switch_to_pe0_level)+1,nys_mg(               &
241                       mg_switch_to_pe0_level)-1:nyn_mg(mg_switch_to_pe0_level)+1,&
242                       nxl_mg(mg_switch_to_pe0_level)-1:nxr_mg(                   &
243                       mg_switch_to_pe0_level)+1) )
244       ELSE
245          ALLOCATE ( p3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
246       ENDIF
247
248       p3 = 0.0_wp
249   
250!
251!--    Ghost boundaries have to be added to divergence array.
252!--    Exchange routine needs to know the grid level!
253       grid_level = maximum_grid_level
254       CALL exchange_horiz( d, 1)
255!
256!--    Set bottom and top boundary conditions
257       d(nzb,:,:) = d(nzb+1,:,:)
258       IF ( ibc_p_t == 1 )  d(nzt+1,:,: ) = d(nzt,:,:)
259!
260!--    Set lateral boundary conditions in non-cyclic case
261       IF ( .NOT. bc_lr_cyc )  THEN
262          IF ( bc_dirichlet_l   .OR.  bc_radiation_l )                         &
263             d(:,:,nxl-1) = d(:,:,nxl)
264          IF ( bc_dirichlet_r   .OR.  bc_radiation_r )                         &
265             d(:,:,nxr+1) = d(:,:,nxr)
266       ENDIF
267       IF ( .NOT. bc_ns_cyc )  THEN
268          IF ( bc_dirichlet_n   .OR.  bc_radiation_n )                         &
269             d(:,nyn+1,:) = d(:,nyn,:)
270          IF ( bc_dirichlet_s   .OR.  bc_radiation_s )                         &
271             d(:,nys-1,:) = d(:,nys,:)
272       ENDIF
273
274!
275!--    Initiation of the multigrid scheme. Does n cycles until the
276!--    residual is smaller than the given limit. The accuracy of the solution
277!--    of the poisson equation will increase with the number of cycles.
278!--    If the number of cycles is preset by the user, this number will be
279!--    carried out regardless of the accuracy.
280       grid_level_count =  0
281       mgcycles         =  0
282       IF ( mg_cycles == -1 )  THEN
283          maximum_mgcycles = 0
284          residual_norm    = 1.0_wp
285       ELSE
286          maximum_mgcycles = mg_cycles
287          residual_norm    = 0.0_wp
288       ENDIF
289
290       DO WHILE ( residual_norm > residual_limit  .OR. &
291                  mgcycles < maximum_mgcycles )
292   
293          CALL next_mg_level_noopt( d, p_loc, p3, r)
294
295!
296!--       Calculate the residual if the user has not preset the number of
297!--       cycles to be performed
298          IF ( maximum_mgcycles == 0 )  THEN
299             CALL resid_noopt( d, p_loc, r )
300             maxerror = SUM( r(nzb+1:nzt,nys:nyn,nxl:nxr)**2 )
301
302#if defined( __parallel )
303             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
304                CALL MPI_ALLREDUCE( maxerror, residual_norm, 1, MPI_REAL,      &
305                                    MPI_SUM, comm2d, ierr)
306#else
307                residual_norm = maxerror
308#endif
309             residual_norm = SQRT( residual_norm )
310          ENDIF
311
312          mgcycles = mgcycles + 1
313
314!
315!--       If the user has not limited the number of cycles, stop the run in case
316!--       of insufficient convergence
317          IF ( mgcycles > 1000  .AND.  mg_cycles == -1 )  THEN
318             message_string = 'no sufficient convergence within 1000 cycles'
319             CALL message( 'poismg_noopt', 'PA0283', 1, 2, 0, 6, 0 )
320          ENDIF
321
322       ENDDO
323
324       DEALLOCATE( p3 )
325
326!
327!--    Unset the grid level. Variable is used to determine the MPI datatypes for
328!--    ghost point exchange
329       grid_level = 0 
330
331       CALL cpu_log( log_point_s(29), 'poismg_noopt', 'stop' )
332
333    END SUBROUTINE poismg_noopt
334
335
336!------------------------------------------------------------------------------!
337! Description:
338! ------------
339!> Computes the residual of the perturbation pressure.
340!------------------------------------------------------------------------------!
341    SUBROUTINE resid_noopt( f_mg, p_mg, r )
342
343
344       USE arrays_3d,                                                          &
345           ONLY:  f1_mg, f2_mg, f3_mg, rho_air_mg
346
347       USE control_parameters,                                                 &
348           ONLY:  bc_lr_cyc, bc_ns_cyc, ibc_p_b, ibc_p_t
349
350       USE grid_variables,                                                     &
351           ONLY:  ddx2_mg, ddy2_mg
352
353       USE indices,                                                            &
354           ONLY:  flags, wall_flags_1, wall_flags_2, wall_flags_3, wall_flags_4,&
355                  wall_flags_5, wall_flags_6, wall_flags_7, wall_flags_8,      &
356                  wall_flags_9, wall_flags_10, nxl_mg, nxr_mg, nys_mg, nyn_mg, &
357                  nzb, nzt_mg
358
359       USE kinds
360
361       IMPLICIT NONE
362
363       INTEGER(iwp) ::  i
364       INTEGER(iwp) ::  j
365       INTEGER(iwp) ::  k
366       INTEGER(iwp) ::  l
367
368       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
369                       nys_mg(grid_level)-1:nyn_mg(grid_level)+1,              &
370                       nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f_mg  !<
371       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
372                       nys_mg(grid_level)-1:nyn_mg(grid_level)+1,              &
373                       nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  p_mg  !<
374       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
375                       nys_mg(grid_level)-1:nyn_mg(grid_level)+1,              &
376                       nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  r     !<
377
378!
379!--    Calculate the residual
380       l = grid_level
381
382!
383!--    Choose flag array of this level
384       SELECT CASE ( l )
385          CASE ( 1 )
386             flags => wall_flags_1
387          CASE ( 2 )
388             flags => wall_flags_2
389          CASE ( 3 )
390             flags => wall_flags_3
391          CASE ( 4 )
392             flags => wall_flags_4
393          CASE ( 5 )
394             flags => wall_flags_5
395          CASE ( 6 )
396             flags => wall_flags_6
397          CASE ( 7 )
398             flags => wall_flags_7
399          CASE ( 8 )
400             flags => wall_flags_8
401          CASE ( 9 )
402             flags => wall_flags_9
403          CASE ( 10 )
404             flags => wall_flags_10
405       END SELECT
406
407!$OMP PARALLEL PRIVATE (i,j,k)
408!$OMP DO
409       DO  i = nxl_mg(l), nxr_mg(l)
410          DO  j = nys_mg(l), nyn_mg(l) 
411             DO  k = nzb+1, nzt_mg(l)
412                r(k,j,i) = f_mg(k,j,i)                                         &
413                           - rho_air_mg(k,l) * ddx2_mg(l) *                    &
414                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
415                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
416                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
417                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
418                           - rho_air_mg(k,l) * ddy2_mg(l) *                    &
419                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
420                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
421                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
422                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
423                           - f2_mg(k,l) *                                      &
424                               ( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
425                                             ( p_mg(k,j,i) - p_mg(k+1,j,i) ) ) &
426                           - f3_mg(k,l) *                                      &
427                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
428                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
429                           + f1_mg(k,l) * p_mg(k,j,i)
430!
431!--             Residual within topography should be zero
432                r(k,j,i) = r(k,j,i) * ( 1.0_wp - IBITS( flags(k,j,i), 6, 1 ) )
433             ENDDO
434          ENDDO
435       ENDDO
436!$OMP END PARALLEL
437
438!
439!--    Horizontal boundary conditions
440       CALL exchange_horiz( r, 1)
441
442       IF ( .NOT. bc_lr_cyc )  THEN
443          IF ( bc_dirichlet_l  .OR.  bc_radiation_l )  THEN
444             r(:,:,nxl_mg(l)-1) = r(:,:,nxl_mg(l))
445          ENDIF
446          IF ( bc_dirichlet_r  .OR.  bc_radiation_r )  THEN
447             r(:,:,nxr_mg(l)+1) = r(:,:,nxr_mg(l))
448          ENDIF
449       ENDIF
450
451       IF ( .NOT. bc_ns_cyc )  THEN
452          IF ( bc_dirichlet_n  .OR.  bc_radiation_n )  THEN
453             r(:,nyn_mg(l)+1,:) = r(:,nyn_mg(l),:)
454          ENDIF
455          IF ( bc_dirichlet_s  .OR.  bc_radiation_s )  THEN
456             r(:,nys_mg(l)-1,:) = r(:,nys_mg(l),:)
457          ENDIF
458       ENDIF
459
460!
461!--    Boundary conditions at bottom and top of the domain.
462!--    These points are not handled by the above loop. Points may be within
463!--    buildings, but that doesn't matter.
464       IF ( ibc_p_b == 1 )  THEN
465          r(nzb,:,: ) = r(nzb+1,:,:)
466       ELSE
467          r(nzb,:,: ) = 0.0_wp
468       ENDIF
469
470       IF ( ibc_p_t == 1 )  THEN
471          r(nzt_mg(l)+1,:,: ) = r(nzt_mg(l),:,:)
472       ELSE
473          r(nzt_mg(l)+1,:,: ) = 0.0_wp
474       ENDIF
475
476
477    END SUBROUTINE resid_noopt
478
479
480!------------------------------------------------------------------------------!
481! Description:
482! ------------
483!> Interpolates the residual on the next coarser grid with "full weighting"
484!> scheme.
485!------------------------------------------------------------------------------!
486    SUBROUTINE restrict_noopt( f_mg, r )
487
488
489       USE control_parameters,                                                 &
490           ONLY:  bc_lr_cyc, bc_ns_cyc, grid_level, ibc_p_b, ibc_p_t
491
492       USE indices,                                                            &
493           ONLY:  flags, wall_flags_1, wall_flags_2, wall_flags_3, wall_flags_4,&
494                  wall_flags_5, wall_flags_6, wall_flags_7, wall_flags_8,      &
495                  wall_flags_9, wall_flags_10, nxl_mg, nxr_mg, nys_mg, nyn_mg, &
496                  nzb, nzt_mg
497
498       USE kinds
499
500       IMPLICIT NONE
501
502       INTEGER(iwp) ::  i    !<
503       INTEGER(iwp) ::  ic   !<
504       INTEGER(iwp) ::  j    !<
505       INTEGER(iwp) ::  jc   !<
506       INTEGER(iwp) ::  k    !<
507       INTEGER(iwp) ::  kc   !<
508       INTEGER(iwp) ::  l    !<
509
510       REAL(wp) ::  rkjim    !<
511       REAL(wp) ::  rkjip    !<
512       REAL(wp) ::  rkjmi    !<
513       REAL(wp) ::  rkjmim   !<
514       REAL(wp) ::  rkjmip   !<
515       REAL(wp) ::  rkjpi    !<
516       REAL(wp) ::  rkjpim   !<
517       REAL(wp) ::  rkjpip   !<
518       REAL(wp) ::  rkmji    !<
519       REAL(wp) ::  rkmjim   !<
520       REAL(wp) ::  rkmjip   !<
521       REAL(wp) ::  rkmjmi   !<
522       REAL(wp) ::  rkmjmim  !<
523       REAL(wp) ::  rkmjmip  !<
524       REAL(wp) ::  rkmjpi   !<
525       REAL(wp) ::  rkmjpim  !<
526       REAL(wp) ::  rkmjpip  !<
527
528       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
529                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
530                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f_mg  !<
531
532       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level+1)+1,                         &
533                           nys_mg(grid_level+1)-1:nyn_mg(grid_level+1)+1,      &
534                           nxl_mg(grid_level+1)-1:nxr_mg(grid_level+1)+1) ::  r !<
535
536!
537!--    Interpolate the residual
538       l = grid_level
539
540!
541!--    Choose flag array of the upper level
542       SELECT CASE ( l+1 )
543          CASE ( 1 )
544             flags => wall_flags_1
545          CASE ( 2 )
546             flags => wall_flags_2
547          CASE ( 3 )
548             flags => wall_flags_3
549          CASE ( 4 )
550             flags => wall_flags_4
551          CASE ( 5 )
552             flags => wall_flags_5
553          CASE ( 6 )
554             flags => wall_flags_6
555          CASE ( 7 )
556             flags => wall_flags_7
557          CASE ( 8 )
558             flags => wall_flags_8
559          CASE ( 9 )
560             flags => wall_flags_9
561          CASE ( 10 )
562             flags => wall_flags_10
563       END SELECT
564       
565!$OMP PARALLEL PRIVATE (i,j,k,ic,jc,kc, rkjim,rkjip,rkjpi,rkjmi,rkjmim,rkjpim, &
566!$OMP rkjmip, rkjpip,rkmji,rkmjim,rkmjip,rkmjpi,rkmjmi,rkmjmim,rkmjpim,rkmjmip,&
567!$OMP rkmjpip          )
568!$OMP DO
569       DO  ic = nxl_mg(l), nxr_mg(l)   
570          i = 2*ic
571          DO  jc = nys_mg(l), nyn_mg(l) 
572             j = 2*jc
573             DO  kc = nzb+1, nzt_mg(l)
574                k = 2*kc-1
575!
576!--             Use implicit Neumann BCs if the respective gridpoint is inside
577!--             the building
578                rkjim   = r(k,j,i-1)     + IBITS( flags(k,j,i-1), 6, 1 ) *     &
579                                           ( r(k,j,i) - r(k,j,i-1) )
580                rkjip   = r(k,j,i+1)     + IBITS( flags(k,j,i+1), 6, 1 ) *     &
581                                           ( r(k,j,i) - r(k,j,i+1) )
582                rkjpi   = r(k,j+1,i)     + IBITS( flags(k,j+1,i), 6, 1 ) *     &
583                                           ( r(k,j,i) - r(k,j+1,i) )
584                rkjmi   = r(k,j-1,i)     + IBITS( flags(k,j-1,i), 6, 1 ) *     &
585                                           ( r(k,j,i) - r(k,j-1,i) )
586                rkjmim  = r(k,j-1,i-1)   + IBITS( flags(k,j-1,i-1), 6, 1 ) *   &
587                                           ( r(k,j,i) - r(k,j-1,i-1) )
588                rkjpim  = r(k,j+1,i-1)   + IBITS( flags(k,j+1,i-1), 6, 1 ) *   &
589                                           ( r(k,j,i) - r(k,j+1,i-1) )
590                rkjmip  = r(k,j-1,i+1)   + IBITS( flags(k,j-1,i+1), 6, 1 ) *   &
591                                           ( r(k,j,i) - r(k,j-1,i+1) )
592                rkjpip  = r(k,j+1,i+1)   + IBITS( flags(k,j+1,i+1), 6, 1 ) *   &
593                                           ( r(k,j,i) - r(k,j+1,i+1) )
594                rkmji   = r(k-1,j,i)     + IBITS( flags(k-1,j,i), 6, 1 ) *     &
595                                           ( r(k,j,i) - r(k-1,j,i) )
596                rkmjim  = r(k-1,j,i-1)   + IBITS( flags(k-1,j,i-1), 6, 1 ) *   &
597                                           ( r(k,j,i) - r(k-1,j,i-1) )
598                rkmjip  = r(k-1,j,i+1)   + IBITS( flags(k-1,j,i+1), 6, 1 ) *   &
599                                           ( r(k,j,i) - r(k-1,j,i+1) )
600                rkmjpi  = r(k-1,j+1,i)   + IBITS( flags(k-1,j+1,i), 6, 1 ) *   &
601                                           ( r(k,j,i) - r(k-1,j+1,i) )
602                rkmjmi  = r(k-1,j-1,i)   + IBITS( flags(k-1,j-1,i), 6, 1 ) *   &
603                                           ( r(k,j,i) - r(k-1,j-1,i) )
604                rkmjmim = r(k-1,j-1,i-1) + IBITS( flags(k-1,j-1,i-1), 6, 1 ) * &
605                                           ( r(k,j,i) - r(k-1,j-1,i-1) )
606                rkmjpim = r(k-1,j+1,i-1) + IBITS( flags(k-1,j+1,i-1), 6, 1 ) * &
607                                           ( r(k,j,i) - r(k-1,j+1,i-1) )
608                rkmjmip = r(k-1,j-1,i+1) + IBITS( flags(k-1,j-1,i+1), 6, 1 ) * &
609                                           ( r(k,j,i) - r(k-1,j-1,i+1) )
610                rkmjpip = r(k-1,j+1,i+1) + IBITS( flags(k-1,j+1,i+1), 6, 1 ) * &
611                                           ( r(k,j,i) - r(k-1,j+1,i+1) )
612
613                f_mg(kc,jc,ic) = 1.0_wp / 64.0_wp * (                         &
614                                 8.0_wp * r(k,j,i)                            &
615                               + 4.0_wp * ( rkjim   + rkjip   +               &
616                                            rkjpi   + rkjmi   )               &
617                               + 2.0_wp * ( rkjmim  + rkjpim  +               &
618                                            rkjmip  + rkjpip  )               &
619                               + 4.0_wp * rkmji                               &
620                               + 2.0_wp * ( rkmjim  + rkmjim  +               &
621                                            rkmjpi  + rkmjmi  )               &
622                               +          ( rkmjmim + rkmjpim +               &
623                                            rkmjmip + rkmjpip )               &
624                               + 4.0_wp * r(k+1,j,i)                          &
625                               + 2.0_wp * ( r(k+1,j,i-1)   + r(k+1,j,i+1)   + &
626                                            r(k+1,j+1,i)   + r(k+1,j-1,i)   ) &
627                               +          ( r(k+1,j-1,i-1) + r(k+1,j+1,i-1) + &
628                                            r(k+1,j-1,i+1) + r(k+1,j+1,i+1) ) &
629                                                    )
630
631!             f_mg(kc,jc,ic) = 1.0_wp / 64.0_wp * (                         &
632!                              8.0_wp * r(k,j,i)                            &
633!                            + 4.0_wp * ( r(k,j,i-1)     + r(k,j,i+1)     + &
634!                                         r(k,j+1,i)     + r(k,j-1,i)     ) &
635!                            + 2.0_wp * ( r(k,j-1,i-1)   + r(k,j+1,i-1)   + &
636!                                         r(k,j-1,i+1)   + r(k,j+1,i+1)   ) &
637!                            + 4.0_wp * r(k-1,j,i)                          &
638!                            + 2.0_wp * ( r(k-1,j,i-1)   + r(k-1,j,i+1)   + &
639!                                         r(k-1,j+1,i)   + r(k-1,j-1,i)   ) &
640!                            +          ( r(k-1,j-1,i-1) + r(k-1,j+1,i-1) + &
641!                                         r(k-1,j-1,i+1) + r(k-1,j+1,i+1) ) &
642!                            + 4.0_wp * r(k+1,j,i)                          &
643!                            + 2.0_wp * ( r(k+1,j,i-1)   + r(k+1,j,i+1)   + &
644!                                         r(k+1,j+1,i)   + r(k+1,j-1,i)   ) &
645!                            +          ( r(k+1,j-1,i-1) + r(k+1,j+1,i-1) + &
646!                                         r(k+1,j-1,i+1) + r(k+1,j+1,i+1) ) &
647!                                                )
648             ENDDO
649          ENDDO
650       ENDDO
651!$OMP END PARALLEL
652
653!
654!--    Horizontal boundary conditions
655       CALL exchange_horiz( f_mg, 1)
656
657       IF ( .NOT. bc_lr_cyc )  THEN
658          IF ( bc_dirichlet_l  .OR.  bc_radiation_l )  THEN
659             f_mg(:,:,nxl_mg(l)-1) = f_mg(:,:,nxl_mg(l))
660          ENDIF
661          IF ( bc_dirichlet_r  .OR.  bc_radiation_r )  THEN
662             f_mg(:,:,nxr_mg(l)+1) = f_mg(:,:,nxr_mg(l))
663          ENDIF
664       ENDIF
665
666       IF ( .NOT. bc_ns_cyc )  THEN
667          IF ( bc_dirichlet_n  .OR.  bc_radiation_n )  THEN
668             f_mg(:,nyn_mg(l)+1,:) = f_mg(:,nyn_mg(l),:)
669          ENDIF
670          IF ( bc_dirichlet_s  .OR.  bc_radiation_s )  THEN
671             f_mg(:,nys_mg(l)-1,:) = f_mg(:,nys_mg(l),:)
672          ENDIF
673       ENDIF
674
675!
676!--    Boundary conditions at bottom and top of the domain.
677!--    These points are not handled by the above loop. Points may be within
678!--    buildings, but that doesn't matter.
679       IF ( ibc_p_b == 1 )  THEN
680          f_mg(nzb,:,: ) = f_mg(nzb+1,:,:)
681       ELSE
682          f_mg(nzb,:,: ) = 0.0_wp
683       ENDIF
684
685       IF ( ibc_p_t == 1 )  THEN
686          f_mg(nzt_mg(l)+1,:,: ) = f_mg(nzt_mg(l),:,:)
687       ELSE
688          f_mg(nzt_mg(l)+1,:,: ) = 0.0_wp
689       ENDIF
690
691
692   END SUBROUTINE restrict_noopt
693
694
695!------------------------------------------------------------------------------!
696! Description:
697! ------------
698!> Interpolates the correction of the perturbation pressure
699!> to the next finer grid.
700!------------------------------------------------------------------------------!
701 SUBROUTINE prolong_noopt( p, temp )
702
703
704    USE control_parameters,                                                    &
705        ONLY:  bc_lr_cyc, bc_ns_cyc, ibc_p_b, ibc_p_t
706    USE indices,                                                               &
707        ONLY:  nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
708
709    USE kinds
710
711    IMPLICIT NONE
712
713    INTEGER(iwp) ::  i  !<
714    INTEGER(iwp) ::  j  !<
715    INTEGER(iwp) ::  k  !<
716    INTEGER(iwp) ::  l  !<
717
718    REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1,                           &
719                        nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1,        &
720                        nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1 ) ::  p  !<
721
722    REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
723                        nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
724                        nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  temp  !<
725
726
727!
728!-- First, store elements of the coarser grid on the next finer grid
729    l = grid_level
730
731!$OMP PARALLEL PRIVATE (i,j,k)
732!$OMP DO
733    DO  i = nxl_mg(l-1), nxr_mg(l-1)
734       DO  j = nys_mg(l-1), nyn_mg(l-1)
735!CDIR NODEP
736          DO  k = nzb+1, nzt_mg(l-1)
737!
738!--          Points of the coarse grid are directly stored on the next finer
739!--          grid
740             temp(2*k-1,2*j,2*i) = p(k,j,i) 
741!
742!--          Points between two coarse-grid points
743             temp(2*k-1,2*j,2*i+1) = 0.5_wp * ( p(k,j,i) + p(k,j,i+1) )
744             temp(2*k-1,2*j+1,2*i) = 0.5_wp * ( p(k,j,i) + p(k,j+1,i) )
745             temp(2*k,2*j,2*i)     = 0.5_wp * ( p(k,j,i) + p(k+1,j,i) )
746!
747!--          Points in the center of the planes stretched by four points
748!--          of the coarse grid cube
749             temp(2*k-1,2*j+1,2*i+1) = 0.25_wp * ( p(k,j,i)   + p(k,j,i+1) + &
750                                                   p(k,j+1,i) + p(k,j+1,i+1) )
751             temp(2*k,2*j,2*i+1)     = 0.25_wp * ( p(k,j,i)   + p(k,j,i+1) + &
752                                                   p(k+1,j,i) + p(k+1,j,i+1) )
753             temp(2*k,2*j+1,2*i)     = 0.25_wp * ( p(k,j,i)   + p(k,j+1,i) + &
754                                                   p(k+1,j,i) + p(k+1,j+1,i) )
755!
756!--          Points in the middle of coarse grid cube
757             temp(2*k,2*j+1,2*i+1) = 0.125_wp * ( p(k,j,i)     + p(k,j,i+1)   + &
758                                                  p(k,j+1,i)   + p(k,j+1,i+1) + &
759                                                  p(k+1,j,i)   + p(k+1,j,i+1) + &
760                                                  p(k+1,j+1,i) + p(k+1,j+1,i+1) )
761          ENDDO
762       ENDDO
763    ENDDO
764!$OMP END PARALLEL
765                         
766!
767!-- Horizontal boundary conditions
768    CALL exchange_horiz( temp, 1)
769
770    IF ( .NOT. bc_lr_cyc )  THEN
771       IF ( bc_dirichlet_l  .OR.  bc_radiation_l )  THEN
772          temp(:,:,nxl_mg(l)-1) = temp(:,:,nxl_mg(l))
773       ENDIF
774       IF ( bc_dirichlet_r  .OR.  bc_radiation_r )  THEN
775          temp(:,:,nxr_mg(l)+1) = temp(:,:,nxr_mg(l))
776       ENDIF
777    ENDIF
778
779    IF ( .NOT. bc_ns_cyc )  THEN
780       IF ( bc_dirichlet_n  .OR.  bc_radiation_n )  THEN
781          temp(:,nyn_mg(l)+1,:) = temp(:,nyn_mg(l),:)
782       ENDIF
783       IF ( bc_dirichlet_s  .OR.  bc_radiation_s )  THEN
784          temp(:,nys_mg(l)-1,:) = temp(:,nys_mg(l),:)
785       ENDIF
786    ENDIF
787
788!
789!-- Bottom and top boundary conditions
790    IF ( ibc_p_b == 1 )  THEN
791       temp(nzb,:,: ) = temp(nzb+1,:,:)
792    ELSE
793       temp(nzb,:,: ) = 0.0_wp
794    ENDIF
795
796    IF ( ibc_p_t == 1 )  THEN
797       temp(nzt_mg(l)+1,:,: ) = temp(nzt_mg(l),:,:)
798    ELSE
799       temp(nzt_mg(l)+1,:,: ) = 0.0_wp
800    ENDIF
801
802 
803 END SUBROUTINE prolong_noopt
804
805
806!------------------------------------------------------------------------------!
807! Description:
808! ------------
809!> Relaxation method for the multigrid scheme. A Gauss-Seidel iteration with
810!> 3D-Red-Black decomposition (GS-RB) is used.
811!------------------------------------------------------------------------------!
812    SUBROUTINE redblack_noopt( f_mg, p_mg )
813
814
815       USE arrays_3d,                                                          &
816           ONLY:  f1_mg, f2_mg, f3_mg, rho_air_mg
817
818       USE control_parameters,                                                 &
819           ONLY:  bc_lr_cyc, bc_ns_cyc, ibc_p_b, ibc_p_t, ngsrb
820
821       USE cpulog,                                                             &
822           ONLY:  cpu_log, log_point_s
823
824       USE grid_variables,                                                     &
825           ONLY:  ddx2_mg, ddy2_mg
826
827       USE indices,                                                            &
828           ONLY:  flags, wall_flags_1, wall_flags_2, wall_flags_3, wall_flags_4,&
829                  wall_flags_5, wall_flags_6, wall_flags_7, wall_flags_8,      &
830                  wall_flags_9, wall_flags_10, nxl_mg, nxr_mg, nys_mg, nyn_mg, &
831                  nzb, nzt_mg
832
833       USE kinds
834
835       IMPLICIT NONE
836
837       INTEGER(iwp) :: color    !<
838       INTEGER(iwp) :: i        !<
839       INTEGER(iwp) :: ic       !<
840       INTEGER(iwp) :: j        !<
841       INTEGER(iwp) :: jc       !<
842       INTEGER(iwp) :: jj       !<
843       INTEGER(iwp) :: k        !<
844       INTEGER(iwp) :: l        !<
845       INTEGER(iwp) :: n        !<
846
847       LOGICAL :: unroll        !<
848
849       REAL(wp) ::  wall_left   !<
850       REAL(wp) ::  wall_north  !<
851       REAL(wp) ::  wall_right  !<
852       REAL(wp) ::  wall_south  !<
853       REAL(wp) ::  wall_total  !<
854       REAL(wp) ::  wall_top    !<
855
856       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
857                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
858                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f_mg  !<
859       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
860                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
861                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  p_mg  !<
862
863       l = grid_level
864
865!
866!--    Choose flag array of this level
867       SELECT CASE ( l )
868          CASE ( 1 )
869             flags => wall_flags_1
870          CASE ( 2 )
871             flags => wall_flags_2
872          CASE ( 3 )
873             flags => wall_flags_3
874          CASE ( 4 )
875             flags => wall_flags_4
876          CASE ( 5 )
877             flags => wall_flags_5
878          CASE ( 6 )
879             flags => wall_flags_6
880          CASE ( 7 )
881             flags => wall_flags_7
882          CASE ( 8 )
883             flags => wall_flags_8
884          CASE ( 9 )
885             flags => wall_flags_9
886          CASE ( 10 )
887             flags => wall_flags_10
888       END SELECT
889
890       unroll = ( MOD( nyn_mg(l)-nys_mg(l)+1, 4 ) == 0  .AND. &
891                  MOD( nxr_mg(l)-nxl_mg(l)+1, 2 ) == 0 )
892
893       DO  n = 1, ngsrb
894         
895          DO  color = 1, 2
896
897             IF ( .NOT. unroll )  THEN
898
899                CALL cpu_log( log_point_s(36), 'redblack_no_unroll_noopt', 'start' )
900
901!
902!--             Without unrolling of loops, no cache optimization
903                DO  i = nxl_mg(l), nxr_mg(l), 2
904                   DO  j = nys_mg(l) + 2 - color, nyn_mg(l), 2 
905                      DO  k = nzb+1, nzt_mg(l), 2
906!                      p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
907!                                ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
908!                              + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
909!                              + f2_mg(k,l) * p_mg(k+1,j,i)                     &
910!                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
911!                                                       )
912
913                         p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                 &
914                             rho_air_mg(k,l) * ddx2_mg(l) *                    &
915                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
916                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
917                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
918                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
919                           + rho_air_mg(k,l) * ddy2_mg(l) *                    &
920                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
921                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
922                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
923                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
924                           + f2_mg(k,l) *                                      &
925                               ( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 7, 1 ) * &
926                                             ( p_mg(k,j,i) - p_mg(k+1,j,i) ) ) &
927                           + f3_mg(k,l) *                                      &
928                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
929                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
930                           - f_mg(k,j,i)                    )
931                      ENDDO
932                   ENDDO
933                ENDDO
934     
935                DO  i = nxl_mg(l)+1, nxr_mg(l), 2
936                   DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
937                      DO  k = nzb+1, nzt_mg(l), 2 
938                         p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                 &
939                               rho_air_mg(k,l) * ddx2_mg(l) *                  &
940                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
941                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
942                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
943                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
944                           + rho_air_mg(k,l) * ddy2_mg(l) *                    &
945                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
946                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
947                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
948                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
949                           + f2_mg(k,l) *                                      &
950                               ( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 7, 1 ) * &
951                                             ( p_mg(k,j,i) - p_mg(k+1,j,i) ) ) &
952                           + f3_mg(k,l) *                                      &
953                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
954                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
955                           - f_mg(k,j,i)                     )
956                      ENDDO
957                   ENDDO
958                ENDDO
959     
960                DO  i = nxl_mg(l), nxr_mg(l), 2
961                   DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
962                      DO  k = nzb+2, nzt_mg(l), 2
963                        p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                  &
964                             rho_air_mg(k,l) * ddx2_mg(l) *                    &
965                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
966                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
967                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
968                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
969                           + rho_air_mg(k,l) * ddy2_mg(l) *                    &
970                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
971                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
972                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
973                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
974                           + f2_mg(k,l) *                                      &
975                               ( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 7, 1 ) * &
976                                             ( p_mg(k,j,i) - p_mg(k+1,j,i) ) ) &
977                           + f3_mg(k,l) *                                      &
978                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
979                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
980                           - f_mg(k,j,i)                    )
981                      ENDDO
982                   ENDDO
983                ENDDO
984
985                DO  i = nxl_mg(l)+1, nxr_mg(l), 2
986                   DO  j = nys_mg(l) + 2 - color, nyn_mg(l), 2
987                      DO  k = nzb+2, nzt_mg(l), 2
988                         p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                 &
989                             rho_air_mg(k,l) * ddx2_mg(l) *                    &
990                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
991                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
992                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
993                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
994                           + rho_air_mg(k,l) * ddy2_mg(l) *                    &
995                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
996                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
997                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
998                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
999                           + f2_mg(k,l) *                                      &
1000                               ( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 7, 1 ) * &
1001                                             ( p_mg(k,j,i) - p_mg(k+1,j,i) ) ) &
1002                           + f3_mg(k,l) *                                      &
1003                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
1004                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
1005                           - f_mg(k,j,i)                    )
1006                      ENDDO
1007                   ENDDO
1008                ENDDO
1009                CALL cpu_log( log_point_s(36), 'redblack_no_unroll_noopt', 'stop' )
1010
1011             ELSE
1012
1013!
1014!--             Loop unrolling along y, only one i loop for better cache use
1015                CALL cpu_log( log_point_s(38), 'redblack_unroll_noopt', 'start' )
1016                DO  ic = nxl_mg(l), nxr_mg(l), 2
1017                   DO  jc = nys_mg(l), nyn_mg(l), 4
1018                      i  = ic
1019                      jj = jc+2-color
1020                      DO  k = nzb+1, nzt_mg(l), 2
1021                         j = jj
1022                         p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                 &
1023                             rho_air_mg(k,l) * ddx2_mg(l) *                    &
1024                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
1025                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
1026                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
1027                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
1028                           + rho_air_mg(k,l) * ddy2_mg(l) *                    &
1029                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
1030                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
1031                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
1032                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
1033                           + f2_mg(k,l) *                                      &
1034                               ( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 7, 1 ) * &
1035                                             ( p_mg(k,j,i) - p_mg(k+1,j,i) ) ) &
1036                           + f3_mg(k,l) *                                      &
1037                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
1038                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
1039                           - f_mg(k,j,i)               )
1040                         j = jj+2
1041                         p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                 &
1042                             rho_air_mg(k,l) * ddx2_mg(l) *                    &
1043                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
1044                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
1045                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
1046                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
1047                           + rho_air_mg(k,l) * ddy2_mg(l) *                    &
1048                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
1049                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
1050                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
1051                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
1052                           + f2_mg(k,l) *                                      &
1053                               ( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 7, 1 ) * &
1054                                             ( p_mg(k,j,i) - p_mg(k+1,j,i) ) ) &
1055                           + f3_mg(k,l) *                                      &
1056                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
1057                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
1058                           - f_mg(k,j,i)                     )
1059                      ENDDO
1060     
1061                      i  = ic+1
1062                      jj = jc+color-1
1063                      DO  k = nzb+1, nzt_mg(l), 2 
1064                         j =jj
1065                         p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                 &
1066                             rho_air_mg(k,l) * ddx2_mg(l) *                    &
1067                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
1068                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
1069                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
1070                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
1071                           + rho_air_mg(k,l) * ddy2_mg(l) *                    &
1072                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
1073                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
1074                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
1075                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
1076                           + f2_mg(k,l) *                                      &
1077                               ( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 7, 1 ) * &
1078                                             ( p_mg(k,j,i) - p_mg(k+1,j,i) ) ) &
1079                           + f3_mg(k,l) *                                      &
1080                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
1081                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
1082                           - f_mg(k,j,i)                     )
1083                         j = jj+2
1084                         p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                &
1085                            rho_air_mg(k,l) * ddx2_mg(l) *                    &
1086                              ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
1087                                            ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
1088                                p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
1089                                            ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
1090                          + rho_air_mg(k,l) * ddy2_mg(l) *                    &
1091                              ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
1092                                            ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
1093                                p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
1094                                            ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
1095                          + f2_mg(k,l) *                                      &
1096                              ( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 7, 1 ) * &
1097                                            ( p_mg(k,j,i) - p_mg(k+1,j,i) ) ) &
1098                          + f3_mg(k,l) *                                      &
1099                              ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
1100                                            ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
1101                          - f_mg(k,j,i)                    )
1102                      ENDDO
1103
1104                      i  = ic
1105                      jj = jc+color-1
1106                      DO  k = nzb+2, nzt_mg(l), 2
1107                         j =jj
1108                         p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                 &
1109                             rho_air_mg(k,l) * ddx2_mg(l) *                    &
1110                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
1111                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
1112                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
1113                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
1114                           + rho_air_mg(k,l) * ddy2_mg(l) *                    &
1115                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
1116                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
1117                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
1118                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
1119                           + f2_mg(k,l) *                                      &
1120                               ( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 7, 1 ) * &
1121                                             ( p_mg(k,j,i) - p_mg(k+1,j,i) ) ) &
1122                           + f3_mg(k,l) *                                      &
1123                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
1124                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
1125                           - f_mg(k,j,i)                    )
1126                         j = jj+2
1127                         p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                 &
1128                             rho_air_mg(k,l) * ddx2_mg(l) *                    &
1129                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
1130                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
1131                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
1132                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
1133                           + rho_air_mg(k,l) * ddy2_mg(l) *                    &
1134                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
1135                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
1136                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
1137                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
1138                           + f2_mg(k,l) *                                      &
1139                               ( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 7, 1 ) * &
1140                                             ( p_mg(k,j,i) - p_mg(k+1,j,i) ) ) &
1141                           + f3_mg(k,l) *                                      &
1142                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
1143                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
1144                           - f_mg(k,j,i)                     )
1145                      ENDDO
1146
1147                      i  = ic+1
1148                      jj = jc+2-color
1149                      DO  k = nzb+2, nzt_mg(l), 2
1150                         j =jj
1151                         p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                 &
1152                             rho_air_mg(k,l) * ddx2_mg(l) *                    &
1153                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
1154                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
1155                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
1156                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
1157                           + rho_air_mg(k,l) * ddy2_mg(l) *                    &
1158                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
1159                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
1160                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
1161                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
1162                           + f2_mg(k,l) *                                      &
1163                               ( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 7, 1 ) * &
1164                                             ( p_mg(k,j,i) - p_mg(k+1,j,i) ) ) &
1165                           + f3_mg(k,l) *                                      &
1166                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
1167                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
1168                           - f_mg(k,j,i)                     )
1169                         j = jj+2
1170                         p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                 &
1171                             rho_air_mg(k,l) * ddx2_mg(l) *                    &
1172                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
1173                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
1174                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
1175                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
1176                           + rho_air_mg(k,l) * ddy2_mg(l) *                    &
1177                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
1178                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
1179                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
1180                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
1181                           + f2_mg(k,l) *                                      &
1182                               ( p_mg(k+1,j,i) + IBITS( flags(k,j,i), 7, 1 ) * &
1183                                             ( p_mg(k,j,i) - p_mg(k+1,j,i) ) ) &
1184                           + f3_mg(k,l) *                                      &
1185                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
1186                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
1187                           - f_mg(k,j,i)                     )
1188                      ENDDO
1189
1190                   ENDDO
1191                ENDDO
1192                CALL cpu_log( log_point_s(38), 'redblack_unroll_noopt', 'stop' )
1193
1194             ENDIF
1195
1196!
1197!--          Horizontal boundary conditions
1198             CALL exchange_horiz( p_mg, 1 )
1199
1200             IF ( .NOT. bc_lr_cyc )  THEN
1201                IF ( bc_dirichlet_l  .OR.  bc_radiation_l )  THEN
1202                   p_mg(:,:,nxl_mg(l)-1) = p_mg(:,:,nxl_mg(l))
1203                ENDIF
1204                IF ( bc_dirichlet_r  .OR.  bc_radiation_r )  THEN
1205                   p_mg(:,:,nxr_mg(l)+1) = p_mg(:,:,nxr_mg(l))
1206                ENDIF
1207             ENDIF
1208
1209             IF ( .NOT. bc_ns_cyc )  THEN
1210                IF ( bc_dirichlet_n  .OR.  bc_radiation_n )  THEN
1211                   p_mg(:,nyn_mg(l)+1,:) = p_mg(:,nyn_mg(l),:)
1212                ENDIF
1213                IF ( bc_dirichlet_s  .OR.  bc_radiation_s )  THEN
1214                   p_mg(:,nys_mg(l)-1,:) = p_mg(:,nys_mg(l),:)
1215                ENDIF
1216             ENDIF
1217
1218!
1219!--          Bottom and top boundary conditions
1220             IF ( ibc_p_b == 1 )  THEN
1221                p_mg(nzb,:,: ) = p_mg(nzb+1,:,:)
1222             ELSE
1223                p_mg(nzb,:,: ) = 0.0_wp
1224             ENDIF
1225
1226             IF ( ibc_p_t == 1 )  THEN
1227                p_mg(nzt_mg(l)+1,:,: ) = p_mg(nzt_mg(l),:,:)
1228             ELSE
1229                p_mg(nzt_mg(l)+1,:,: ) = 0.0_wp
1230             ENDIF
1231
1232          ENDDO
1233
1234       ENDDO
1235
1236!
1237!--    Set pressure within topography and at the topography surfaces
1238!$OMP PARALLEL PRIVATE (i,j,k,wall_left,wall_north,wall_right,wall_south,wall_top,wall_total)
1239!$OMP DO
1240       DO  i = nxl_mg(l), nxr_mg(l)
1241          DO  j = nys_mg(l), nyn_mg(l) 
1242             DO  k = nzb, nzt_mg(l)
1243!
1244!--             First, set pressure inside topography to zero
1245                p_mg(k,j,i) = p_mg(k,j,i) * ( 1.0_wp - IBITS( flags(k,j,i), 6, 1 ) )
1246!
1247!--             Second, determine if the gridpoint inside topography is adjacent
1248!--             to a wall and set its value to a value given by the average of
1249!--             those values obtained from Neumann boundary condition
1250                wall_left  = IBITS( flags(k,j,i-1), 5, 1 )
1251                wall_right = IBITS( flags(k,j,i+1), 4, 1 )
1252                wall_south = IBITS( flags(k,j-1,i), 3, 1 )
1253                wall_north = IBITS( flags(k,j+1,i), 2, 1 )
1254                wall_top   = IBITS( flags(k+1,j,i), 0, 1 )
1255                wall_total = wall_left + wall_right + wall_south + wall_north + &
1256                             wall_top
1257
1258                IF ( wall_total > 0.0_wp )  THEN
1259                   p_mg(k,j,i) = 1.0_wp / wall_total *                 &
1260                                        ( wall_left  * p_mg(k,j,i-1) + &
1261                                          wall_right * p_mg(k,j,i+1) + &
1262                                          wall_south * p_mg(k,j-1,i) + &
1263                                          wall_north * p_mg(k,j+1,i) + &
1264                                          wall_top   * p_mg(k+1,j,i) )
1265                ENDIF
1266             ENDDO
1267          ENDDO
1268       ENDDO
1269!$OMP END PARALLEL
1270
1271!
1272!--    One more time horizontal boundary conditions
1273       CALL exchange_horiz( p_mg, 1)
1274
1275
1276    END SUBROUTINE redblack_noopt
1277
1278
1279
1280!------------------------------------------------------------------------------!
1281! Description:
1282! ------------
1283!> Gather subdomain data from all PEs.
1284!------------------------------------------------------------------------------!
1285    SUBROUTINE mg_gather_noopt( f2, f2_sub )
1286
1287       USE cpulog,                                                             &
1288           ONLY:  cpu_log, log_point_s
1289
1290       USE indices,                                                            &
1291           ONLY:  mg_loc_ind, nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
1292
1293       USE kinds
1294
1295       USE pegrid
1296
1297       IMPLICIT NONE
1298
1299       INTEGER(iwp) ::  i       !<
1300       INTEGER(iwp) ::  il      !<
1301       INTEGER(iwp) ::  ir      !<
1302       INTEGER(iwp) ::  j       !<
1303       INTEGER(iwp) ::  jn      !<
1304       INTEGER(iwp) ::  js      !<
1305       INTEGER(iwp) ::  k       !<
1306       INTEGER(iwp) ::  nwords  !<
1307
1308       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
1309                       nys_mg(grid_level)-1:nyn_mg(grid_level)+1,              &
1310                       nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f2    !<
1311       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
1312                       nys_mg(grid_level)-1:nyn_mg(grid_level)+1,              &
1313                       nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f2_l  !<
1314
1315       REAL(wp), DIMENSION(nzb:mg_loc_ind(5,myid)+1,                           &
1316                           mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1,          &
1317                           mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) ::  f2_sub  !<
1318
1319
1320#if defined( __parallel )
1321       CALL cpu_log( log_point_s(34), 'mg_gather_noopt', 'start' )
1322
1323       f2_l = 0.0_wp
1324
1325!
1326!--    Store the local subdomain array on the total array
1327       js = mg_loc_ind(3,myid)
1328       IF ( south_border_pe )  js = js - 1
1329       jn = mg_loc_ind(4,myid)
1330       IF ( north_border_pe )  jn = jn + 1
1331       il = mg_loc_ind(1,myid)
1332       IF ( left_border_pe )   il = il - 1
1333       ir = mg_loc_ind(2,myid)
1334       IF ( right_border_pe )  ir = ir + 1
1335       DO  i = il, ir
1336          DO  j = js, jn
1337             DO  k = nzb, nzt_mg(grid_level)+1
1338                f2_l(k,j,i) = f2_sub(k,j,i)
1339             ENDDO
1340          ENDDO
1341       ENDDO
1342
1343!
1344!--    Find out the number of array elements of the total array
1345       nwords = SIZE( f2 )
1346
1347!
1348!--    Gather subdomain data from all PEs
1349       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1350       CALL MPI_ALLREDUCE( f2_l(nzb,nys_mg(grid_level)-1,nxl_mg(grid_level)-1),&
1351                           f2(nzb,nys_mg(grid_level)-1,nxl_mg(grid_level)-1),  &
1352                           nwords, MPI_REAL, MPI_SUM, comm2d, ierr )
1353
1354       CALL cpu_log( log_point_s(34), 'mg_gather_noopt', 'stop' )
1355#endif
1356       
1357    END SUBROUTINE mg_gather_noopt
1358
1359
1360
1361!------------------------------------------------------------------------------!
1362! Description:
1363! ------------
1364!> @todo It may be possible to improve the speed of this routine by using
1365!>       non-blocking communication
1366!------------------------------------------------------------------------------!
1367    SUBROUTINE mg_scatter_noopt( p2, p2_sub )
1368
1369       USE cpulog,                                                             &
1370           ONLY:  cpu_log, log_point_s
1371
1372       USE indices,                                                            &
1373           ONLY:  mg_loc_ind, nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
1374
1375       USE kinds
1376
1377       USE pegrid
1378
1379       IMPLICIT NONE
1380
1381       INTEGER(iwp) ::  nwords  !<
1382
1383       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1,                         &
1384                           nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1,      &
1385                           nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) ::  p2  !<
1386
1387       REAL(wp), DIMENSION(nzb:mg_loc_ind(5,myid)+1,                           &
1388                           mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1,          &
1389                           mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) ::  p2_sub  !<
1390
1391!
1392!--    Find out the number of array elements of the subdomain array
1393       nwords = SIZE( p2_sub )
1394
1395#if defined( __parallel )
1396       CALL cpu_log( log_point_s(35), 'mg_scatter_noopt', 'start' )
1397
1398       p2_sub = p2(:,mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, &
1399                     mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1)
1400
1401       CALL cpu_log( log_point_s(35), 'mg_scatter_noopt', 'stop' )
1402#endif
1403       
1404    END SUBROUTINE mg_scatter_noopt
1405
1406
1407!------------------------------------------------------------------------------!
1408! Description:
1409! ------------
1410!> This is where the multigrid technique takes place. V- and W- Cycle are
1411!> implemented and steered by the parameter "gamma". Parameter "nue" determines
1412!> the convergence of the multigrid iterative solution. There are nue times
1413!> RB-GS iterations. It should be set to "1" or "2", considering the time effort
1414!> one would like to invest. Last choice shows a very good converging factor,
1415!> but leads to an increase in computing time.
1416!------------------------------------------------------------------------------!
1417    RECURSIVE SUBROUTINE next_mg_level_noopt( f_mg, p_mg, p3, r )
1418
1419       USE control_parameters,                                                 &
1420           ONLY:  bc_lr_dirrad, bc_lr_raddir, bc_ns_dirrad, bc_ns_raddir,      &
1421                  gamma_mg, grid_level_count, ibc_p_b, ibc_p_t,                &
1422                  maximum_grid_level,                                          &
1423                  mg_switch_to_pe0_level, mg_switch_to_pe0, ngsrb
1424
1425
1426       USE indices,                                                            &
1427           ONLY:  mg_loc_ind, nxl, nxl_mg, nxr, nxr_mg, nys, nys_mg, nyn,      &
1428                  nyn_mg, nzb, nzt, nzt_mg
1429
1430       USE kinds
1431
1432       USE pegrid
1433
1434       IMPLICIT NONE
1435
1436       INTEGER(iwp) ::  i            !<
1437       INTEGER(iwp) ::  j            !<
1438       INTEGER(iwp) ::  k            !<
1439       INTEGER(iwp) ::  nxl_mg_save  !<
1440       INTEGER(iwp) ::  nxr_mg_save  !<
1441       INTEGER(iwp) ::  nyn_mg_save  !<
1442       INTEGER(iwp) ::  nys_mg_save  !<
1443       INTEGER(iwp) ::  nzt_mg_save  !<
1444
1445       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
1446                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
1447                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg  !<
1448       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
1449                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
1450                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p_mg  !<
1451       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
1452                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
1453                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p3    !<
1454       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
1455                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
1456                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: r     !<
1457
1458       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1,                         &
1459                           nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1,      &
1460                           nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) ::  f2  !<
1461       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1,                         &
1462                           nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1,      &
1463                           nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) ::  p2  !<
1464
1465       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  f2_sub  !<
1466       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  p2_sub  !<
1467
1468!
1469!--    Restriction to the coarsest grid
1470    10 IF ( grid_level == 1 )  THEN
1471
1472!
1473!--       Solution on the coarsest grid. Double the number of Gauss-Seidel
1474!--       iterations in order to get a more accurate solution.
1475          ngsrb = 2 * ngsrb
1476
1477          CALL redblack_noopt( f_mg, p_mg )
1478
1479          ngsrb = ngsrb / 2
1480
1481
1482       ELSEIF ( grid_level /= 1 )  THEN
1483
1484          grid_level_count(grid_level) = grid_level_count(grid_level) + 1
1485
1486!
1487!--       Solution on the actual grid level
1488          CALL redblack_noopt( f_mg, p_mg )
1489
1490!
1491!--       Determination of the actual residual
1492          CALL resid_noopt( f_mg, p_mg, r )
1493
1494!
1495!--       Restriction of the residual (finer grid values!) to the next coarser
1496!--       grid. Therefore, the grid level has to be decremented now. nxl..nzt have
1497!--       to be set to the coarse grid values, because these variables are needed
1498!--       for the exchange of ghost points in routine exchange_horiz
1499          grid_level = grid_level - 1
1500          nxl = nxl_mg(grid_level)
1501          nys = nys_mg(grid_level)
1502          nxr = nxr_mg(grid_level)
1503          nyn = nyn_mg(grid_level)
1504          nzt = nzt_mg(grid_level)
1505
1506          IF ( grid_level == mg_switch_to_pe0_level )  THEN
1507
1508!
1509!--          From this level on, calculations are done on PE0 only.
1510!--          First, carry out restriction on the subdomain.
1511!--          Therefore, indices of the level have to be changed to subdomain values
1512!--          in between (otherwise, the restrict routine would expect
1513!--          the gathered array)
1514
1515             nxl_mg_save = nxl_mg(grid_level)
1516             nxr_mg_save = nxr_mg(grid_level)
1517             nys_mg_save = nys_mg(grid_level)
1518             nyn_mg_save = nyn_mg(grid_level)
1519             nzt_mg_save = nzt_mg(grid_level)
1520             nxl_mg(grid_level) = mg_loc_ind(1,myid)
1521             nxr_mg(grid_level) = mg_loc_ind(2,myid)
1522             nys_mg(grid_level) = mg_loc_ind(3,myid)
1523             nyn_mg(grid_level) = mg_loc_ind(4,myid)
1524             nzt_mg(grid_level) = mg_loc_ind(5,myid)
1525             nxl = mg_loc_ind(1,myid)
1526             nxr = mg_loc_ind(2,myid)
1527             nys = mg_loc_ind(3,myid)
1528             nyn = mg_loc_ind(4,myid)
1529             nzt = mg_loc_ind(5,myid)
1530
1531             ALLOCATE( f2_sub(nzb:nzt_mg(grid_level)+1,                    &
1532                              nys_mg(grid_level)-1:nyn_mg(grid_level)+1,   &
1533                              nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) )
1534
1535             CALL restrict_noopt( f2_sub, r )
1536
1537!
1538!--          Restore the correct indices of this level
1539             nxl_mg(grid_level) = nxl_mg_save
1540             nxr_mg(grid_level) = nxr_mg_save
1541             nys_mg(grid_level) = nys_mg_save
1542             nyn_mg(grid_level) = nyn_mg_save
1543             nzt_mg(grid_level) = nzt_mg_save
1544             nxl = nxl_mg(grid_level)
1545             nxr = nxr_mg(grid_level)
1546             nys = nys_mg(grid_level)
1547             nyn = nyn_mg(grid_level)
1548             nzt = nzt_mg(grid_level)
1549!
1550!--          Gather all arrays from the subdomains on PE0
1551             CALL mg_gather_noopt( f2, f2_sub )
1552
1553!
1554!--          Set switch for routine exchange_horiz, that no ghostpoint exchange
1555!--          has to be carried out from now on
1556             mg_switch_to_pe0 = .TRUE.
1557
1558!
1559!--          In case of non-cyclic lateral boundary conditions, both in- and
1560!--          outflow conditions have to be used on all PEs after the switch,
1561!--          because then they have the total domain.
1562             IF ( bc_lr_dirrad )  THEN
1563                bc_dirichlet_l  = .TRUE.
1564                bc_dirichlet_r  = .FALSE.
1565                bc_radiation_l = .FALSE.
1566                bc_radiation_r = .TRUE.
1567             ELSEIF ( bc_lr_raddir )  THEN
1568                bc_dirichlet_l  = .FALSE.
1569                bc_dirichlet_r  = .TRUE.
1570                bc_radiation_l = .TRUE.
1571                bc_radiation_r = .FALSE.
1572             ELSEIF ( child_domain  .OR.  nesting_offline )  THEN
1573                bc_dirichlet_l = .TRUE.
1574                bc_dirichlet_r = .TRUE.
1575             ENDIF
1576
1577             IF ( bc_ns_dirrad )  THEN
1578                bc_dirichlet_n  = .TRUE.
1579                bc_dirichlet_s  = .FALSE.
1580                bc_radiation_n = .FALSE.
1581                bc_radiation_s = .TRUE.
1582             ELSEIF ( bc_ns_raddir )  THEN
1583                bc_dirichlet_n  = .FALSE.
1584                bc_dirichlet_s  = .TRUE.
1585                bc_radiation_n = .TRUE.
1586                bc_radiation_s = .FALSE.
1587             ELSEIF ( child_domain  .OR.  nesting_offline )  THEN
1588                bc_dirichlet_s = .TRUE.
1589                bc_dirichlet_n = .TRUE.
1590             ENDIF
1591
1592             DEALLOCATE( f2_sub )
1593
1594          ELSE
1595
1596             CALL restrict_noopt( f2, r )
1597
1598          ENDIF
1599
1600          p2 = 0.0_wp
1601
1602!
1603!--       Repeat the same procedure till the coarsest grid is reached
1604          CALL next_mg_level_noopt( f2, p2, p3, r )
1605
1606       ENDIF
1607
1608!
1609!--    Now follows the prolongation
1610       IF ( grid_level >= 2 )  THEN
1611
1612!
1613!--       Prolongation of the new residual. The values are transferred
1614!--       from the coarse to the next finer grid.
1615          IF ( grid_level == mg_switch_to_pe0_level+1 )  THEN
1616
1617#if defined( __parallel )
1618!
1619!--          At this level, the new residual first has to be scattered from
1620!--          PE0 to the other PEs
1621             ALLOCATE( p2_sub(nzb:mg_loc_ind(5,myid)+1,             &
1622                       mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1,   &
1623                       mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) )
1624
1625             CALL mg_scatter_noopt( p2, p2_sub )
1626
1627!
1628!--          Therefore, indices of the previous level have to be changed to
1629!--          subdomain values in between (otherwise, the prolong routine would
1630!--          expect the gathered array)
1631             nxl_mg_save = nxl_mg(grid_level-1)
1632             nxr_mg_save = nxr_mg(grid_level-1)
1633             nys_mg_save = nys_mg(grid_level-1)
1634             nyn_mg_save = nyn_mg(grid_level-1)
1635             nzt_mg_save = nzt_mg(grid_level-1)
1636             nxl_mg(grid_level-1) = mg_loc_ind(1,myid)
1637             nxr_mg(grid_level-1) = mg_loc_ind(2,myid)
1638             nys_mg(grid_level-1) = mg_loc_ind(3,myid)
1639             nyn_mg(grid_level-1) = mg_loc_ind(4,myid)
1640             nzt_mg(grid_level-1) = mg_loc_ind(5,myid)
1641
1642!
1643!--          Set switch for routine exchange_horiz, that ghostpoint exchange
1644!--          has to be carried again out from now on
1645             mg_switch_to_pe0 = .FALSE.
1646
1647!
1648!--          For non-cyclic lateral boundary conditions and in case of nesting,
1649!--          restore the in-/outflow conditions.
1650             bc_dirichlet_l = .FALSE.;  bc_dirichlet_r = .FALSE.
1651             bc_dirichlet_n = .FALSE.;  bc_dirichlet_s = .FALSE.
1652             bc_radiation_l = .FALSE.;  bc_radiation_r = .FALSE.
1653             bc_radiation_n = .FALSE.;  bc_radiation_s = .FALSE.
1654
1655             IF ( pleft == MPI_PROC_NULL )  THEN
1656                IF ( bc_lr_dirrad  .OR.  child_domain  .OR.  nesting_offline )  &
1657                THEN
1658                   bc_dirichlet_l = .TRUE.
1659                ELSEIF ( bc_lr_raddir )  THEN
1660                   bc_radiation_l = .TRUE.
1661                ENDIF
1662             ENDIF
1663
1664             IF ( pright == MPI_PROC_NULL )  THEN
1665                IF ( bc_lr_dirrad )  THEN
1666                   bc_radiation_r = .TRUE.
1667                ELSEIF ( bc_lr_raddir  .OR.  child_domain  .OR.                 &
1668                         nesting_offline )  THEN
1669                   bc_dirichlet_r = .TRUE.
1670                ENDIF
1671             ENDIF
1672
1673             IF ( psouth == MPI_PROC_NULL )  THEN
1674                IF ( bc_ns_dirrad )  THEN
1675                   bc_radiation_s = .TRUE.
1676                ELSEIF ( bc_ns_raddir  .OR.  child_domain  .OR.                 &
1677                         nesting_offline )  THEN
1678                   bc_dirichlet_s = .TRUE.
1679                ENDIF
1680             ENDIF
1681
1682             IF ( pnorth == MPI_PROC_NULL )  THEN
1683                IF ( bc_ns_dirrad  .OR.  child_domain  .OR.  nesting_offline )  &
1684                THEN
1685                   bc_dirichlet_n = .TRUE.
1686                ELSEIF ( bc_ns_raddir )  THEN
1687                   bc_radiation_n = .TRUE.
1688                ENDIF
1689             ENDIF
1690
1691             CALL prolong_noopt( p2_sub, p3 )
1692
1693!
1694!--          Restore the correct indices of the previous level
1695             nxl_mg(grid_level-1) = nxl_mg_save
1696             nxr_mg(grid_level-1) = nxr_mg_save
1697             nys_mg(grid_level-1) = nys_mg_save
1698             nyn_mg(grid_level-1) = nyn_mg_save
1699             nzt_mg(grid_level-1) = nzt_mg_save
1700
1701             DEALLOCATE( p2_sub )
1702#endif
1703
1704          ELSE
1705
1706             CALL prolong_noopt( p2, p3 )
1707
1708          ENDIF
1709
1710!
1711!--       Computation of the new pressure correction. Therefore,
1712!--       values from prior grids are added up automatically stage by stage.
1713          DO  i = nxl_mg(grid_level)-1, nxr_mg(grid_level)+1
1714             DO  j = nys_mg(grid_level)-1, nyn_mg(grid_level)+1
1715                DO  k = nzb, nzt_mg(grid_level)+1 
1716                   p_mg(k,j,i) = p_mg(k,j,i) + p3(k,j,i)
1717                ENDDO
1718             ENDDO
1719          ENDDO
1720
1721!
1722!--       Relaxation of the new solution
1723          CALL redblack_noopt( f_mg, p_mg )
1724
1725       ENDIF
1726
1727
1728!
1729!--    The following few lines serve the steering of the multigrid scheme
1730       IF ( grid_level == maximum_grid_level )  THEN
1731
1732          GOTO 20
1733
1734       ELSEIF ( grid_level /= maximum_grid_level  .AND.  grid_level /= 1  .AND. &
1735                grid_level_count(grid_level) /= gamma_mg )  THEN
1736
1737          GOTO 10
1738
1739       ENDIF
1740
1741!
1742!--    Reset counter for the next call of poismg_noopt
1743       grid_level_count(grid_level) = 0
1744
1745!
1746!--    Continue with the next finer level. nxl..nzt have to be
1747!--    set to the finer grid values, because these variables are needed for the
1748!--    exchange of ghost points in routine exchange_horiz
1749       grid_level = grid_level + 1
1750       nxl = nxl_mg(grid_level)
1751       nxr = nxr_mg(grid_level)
1752       nys = nys_mg(grid_level)
1753       nyn = nyn_mg(grid_level)
1754       nzt = nzt_mg(grid_level)
1755
1756    20 CONTINUE
1757
1758    END SUBROUTINE next_mg_level_noopt
1759
1760
1761
1762    SUBROUTINE poismg_noopt_init
1763
1764       USE control_parameters,                                                 &
1765           ONLY:  bc_lr_cyc, bc_ns_cyc, masking_method, maximum_grid_level
1766
1767       USE indices,                                                            &
1768           ONLY:  flags, nxl_mg, nxr_mg, nyn_mg, nys_mg, nzb, nzt_mg,          &
1769                  wall_flags_0, wall_flags_1,                                  &
1770                  wall_flags_10, wall_flags_2, wall_flags_3,  wall_flags_4,    &
1771                  wall_flags_5, wall_flags_6, wall_flags_7, wall_flags_8,      &
1772                  wall_flags_9
1773
1774       IMPLICIT NONE
1775
1776       INTEGER(iwp) ::  i             !< index variable along x
1777       INTEGER(iwp) ::  inc           !< incremental parameter for coarsening grid level
1778       INTEGER(iwp) ::  j             !< index variable along y
1779       INTEGER(iwp) ::  k             !< index variable along z
1780       INTEGER(iwp) ::  l             !< loop variable indication current grid level
1781       INTEGER(iwp) ::  nxl_l         !< index of left PE boundary for multigrid level
1782       INTEGER(iwp) ::  nxr_l         !< index of right PE boundary for multigrid level
1783       INTEGER(iwp) ::  nyn_l         !< index of north PE boundary for multigrid level
1784       INTEGER(iwp) ::  nys_l         !< index of south PE boundary for multigrid level
1785       INTEGER(iwp) ::  nzt_l         !< index of top PE boundary for multigrid level
1786
1787       INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE ::  topo_tmp
1788!
1789!--    Gridpoint increment of the current level.
1790       inc = 1
1791       DO  l = maximum_grid_level, 1 , -1
1792!
1793!--       Set grid_level as it is required for exchange_horiz_2d_int
1794          grid_level = l
1795
1796          nxl_l = nxl_mg(l)
1797          nxr_l = nxr_mg(l)
1798          nys_l = nys_mg(l)
1799          nyn_l = nyn_mg(l)
1800          nzt_l = nzt_mg(l)
1801!
1802!--       Assign the flag level to be calculated
1803          SELECT CASE ( l )
1804             CASE ( 1 )
1805                flags => wall_flags_1
1806             CASE ( 2 )
1807                flags => wall_flags_2
1808             CASE ( 3 )
1809                flags => wall_flags_3
1810             CASE ( 4 )
1811                flags => wall_flags_4
1812             CASE ( 5 )
1813                flags => wall_flags_5
1814             CASE ( 6 )
1815                flags => wall_flags_6
1816             CASE ( 7 )
1817                flags => wall_flags_7
1818             CASE ( 8 )
1819                flags => wall_flags_8
1820             CASE ( 9 )
1821                flags => wall_flags_9
1822             CASE ( 10 )
1823                flags => wall_flags_10
1824          END SELECT
1825
1826!
1827!--       Depending on the grid level, set the respective bits in case of
1828!--       neighbouring walls
1829!--       Bit 0:  wall to the bottom
1830!--       Bit 1:  wall to the top (not realized in remaining PALM code so far)
1831!--       Bit 2:  wall to the south
1832!--       Bit 3:  wall to the north
1833!--       Bit 4:  wall to the left
1834!--       Bit 5:  wall to the right
1835!--       Bit 6:  inside building
1836
1837          flags = 0
1838!
1839!--       In case of masking method, flags are not set and multigrid method
1840!--       works like FFT-solver
1841          IF ( .NOT. masking_method )  THEN
1842
1843!
1844!--          Allocate temporary array for topography heights on coarser grid
1845!--          level. Please note, 2 ghoist points are required, in order to
1846!--          calculate flags() on the interior ghost point.
1847             ALLOCATE( topo_tmp(nzb:nzt_l+1,nys_l-1:nyn_l+1,nxl_l-1:nxr_l+1) )
1848             topo_tmp = 0
1849             
1850             DO  i = nxl_l, nxr_l
1851                DO  j = nys_l, nyn_l
1852                   DO  k = nzb, nzt_l
1853                      topo_tmp(k,j,i) = wall_flags_0(k*inc,j*inc,i*inc)
1854                   ENDDO
1855                ENDDO
1856             ENDDO
1857             topo_tmp(nzt_l+1,:,:) = topo_tmp(nzt_l,:,:)
1858!
1859!--          Exchange ghost points on respective multigrid level. 2 ghost points
1860!--          are required, in order to calculate flags on
1861!--          nys_l-1 / nyn_l+1 / nxl_l-1 / nxr_l+1.
1862             CALL exchange_horiz_int( topo_tmp, nys_l, nyn_l, nxl_l, nxr_l, nzt_l, 1 )
1863!
1864!--          Set non-cyclic boundary conditions on respective multigrid level
1865             IF ( .NOT. bc_ns_cyc )  THEN
1866                IF ( bc_dirichlet_s  .OR.  bc_radiation_s )  THEN
1867!                    topo_tmp(:,-2,:) = topo_tmp(:,0,:)
1868                   topo_tmp(:,-1,:) = topo_tmp(:,0,:)
1869                ENDIF
1870                IF ( bc_dirichlet_n  .OR.  bc_radiation_n )  THEN
1871!                    topo_tmp(:,nyn_l+2,:) = topo_tmp(:,nyn_l,:)
1872                   topo_tmp(:,nyn_l+1,:) = topo_tmp(:,nyn_l,:)
1873                ENDIF
1874             ENDIF
1875             IF ( .NOT. bc_lr_cyc )  THEN
1876                IF ( bc_dirichlet_l  .OR.  bc_radiation_l )  THEN
1877!                    topo_tmp(:,:,-2) = topo_tmp(:,:,0)
1878                   topo_tmp(:,:,-1) = topo_tmp(:,:,0)
1879                ENDIF
1880                IF ( bc_dirichlet_r  .OR.  bc_radiation_r )  THEN
1881!                    topo_tmp(:,:,nxr_l+2) = topo_tmp(:,:,nxr_l)     
1882                   topo_tmp(:,:,nxr_l+1) = topo_tmp(:,:,nxr_l)   
1883                ENDIF       
1884             ENDIF
1885                       
1886             DO  i = nxl_l, nxr_l
1887                DO  j = nys_l, nyn_l
1888                   DO  k = nzb, nzt_l     
1889!
1890!--                   Inside/outside building (inside building does not need
1891!--                   further tests for walls)
1892                      IF ( .NOT. BTEST( topo_tmp(k,j,i), 0 ) )  THEN
1893
1894                         flags(k,j,i) = IBSET( flags(k,j,i), 6 )
1895
1896                      ELSE
1897!
1898!--                      Bottom wall
1899                         IF ( .NOT. BTEST( topo_tmp(k-1,j,i), 0 ) )  THEN
1900                            flags(k,j,i) = IBSET( flags(k,j,i), 0 )
1901                         ENDIF
1902!
1903!--                      South wall
1904                         IF ( .NOT. BTEST( topo_tmp(k,j-1,i), 0 ) )  THEN
1905                            flags(k,j,i) = IBSET( flags(k,j,i), 2 )
1906                         ENDIF
1907!
1908!--                      North wall
1909                         IF ( .NOT. BTEST( topo_tmp(k,j+1,i), 0 ) )  THEN
1910                            flags(k,j,i) = IBSET( flags(k,j,i), 3 )
1911                         ENDIF
1912!
1913!--                      Left wall
1914                         IF ( .NOT. BTEST( topo_tmp(k,j,i-1), 0 ) )  THEN
1915                            flags(k,j,i) = IBSET( flags(k,j,i), 4 )
1916                         ENDIF
1917!
1918!--                      Right wall
1919                         IF ( .NOT. BTEST( topo_tmp(k,j,i+1), 0 ) )  THEN
1920                            flags(k,j,i) = IBSET( flags(k,j,i), 5 )
1921                         ENDIF
1922!
1923!--                      Top wall
1924                         IF ( .NOT. BTEST( topo_tmp(k+1,j,i), 0 ) )  THEN
1925                            flags(k,j,i) = IBSET( flags(k,j,i), 7 )
1926                         ENDIF
1927
1928                      ENDIF
1929                           
1930                   ENDDO
1931                ENDDO
1932             ENDDO
1933             flags(nzt_l+1,:,:) = flags(nzt_l,:,:)
1934
1935             CALL exchange_horiz_int( flags, nys_l, nyn_l, nxl_l, nxr_l, nzt_l, 1 )
1936
1937             DEALLOCATE( topo_tmp )
1938
1939          ENDIF
1940
1941          inc = inc * 2
1942
1943       ENDDO
1944!
1945!--    Reset grid_level to "normal" grid
1946       grid_level = 0
1947
1948    END SUBROUTINE poismg_noopt_init
1949
1950 END MODULE poismg_noopt_mod
Note: See TracBrowser for help on using the repository browser.