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

Last change on this file since 4023 was 3655, checked in by knoop, 5 years ago

Bugfix: made "unit" and "found" intend INOUT in module interface subroutines + automatic copyright update

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