source: palm/trunk/SOURCE/poismg_noopt.f90 @ 2037

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

Anelastic approximation implemented

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