source: palm/trunk/SOURCE/poismg.f90 @ 1

Last change on this file since 1 was 1, checked in by raasch, 15 years ago

Initial repository layout and content

File size: 39.9 KB
Line 
1 SUBROUTINE poismg( r )
2
3!------------------------------------------------------------------------------!
4! Attention: Loop unrolling and cache optimization in SOR-Red/Black method
5!            still does not bring the expected speedup on ibm! Further work
6!            is required.
7!
8! Actual revisions:
9! -----------------
10!
11!
12! Former revisions:
13! -----------------
14! $Log: poismg.f90,v $
15! Revision 1.6  2005/03/26 20:55:54  raasch
16! Implementation of non-cyclic (Neumann) horizontal boundary conditions,
17! routine prolong simplified (one call of exchange_horiz spared)
18!
19! Revision 1.5  2003/10/29 09:04:34  raasch
20! On a defined level, data are gathered on PE0 and further calculations are
21! carried out only on this PE.
22!
23! Revision 1.4  2003/03/16 09:45:57  raasch
24! Two underscores (_) are placed in front of all define-strings
25!
26! Revision 1.3  2001/07/23 12:00:31  raasch
27! Neumann boundary conditions for the upper boundary and Diriclet conditions
28! for the lower boundary implemented
29!
30! Revision 1.2  2001/07/20 15:21:55  raasch
31! Praeprocessor directives included for allreduce
32!
33! Revision 1.1  2001/07/20 13:10:51  raasch
34! Initial revision
35!
36!
37! Description:
38! ------------
39! Solves the Poisson equation for the perturbation pressure with a multigrid
40! V- or W-Cycle scheme.
41!
42! This multigrid method was originally developed for PALM by Joerg Uhlenbrock,
43! September 2000 - July 2001.
44!------------------------------------------------------------------------------!
45
46    USE arrays_3d
47    USE control_parameters
48    USE cpulog   
49    USE grid_variables
50    USE indices
51    USE interfaces
52    USE pegrid
53
54    IMPLICIT NONE
55
56    REAL    ::  maxerror, maximum_mgcycles, residual_norm
57
58    REAL, DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) ::  r
59
60    REAL, DIMENSION(:,:,:), ALLOCATABLE ::  p3
61
62
63    CALL cpu_log( log_point_s(29), 'poismg', 'start' )
64
65
66!
67!-- Initialize arrays and variables used in this subroutine
68    ALLOCATE ( p3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
69
70
71!
72!-- Some boundaries have to be added to divergence array
73    CALL exchange_horiz( d, 0, 0 )
74    d(nzb,:,:) = d(nzb+1,:,:)
75
76!
77!-- Initiation of the multigrid scheme. Does n cycles until the
78!-- residual is smaller than the given limit. The accuracy of the solution
79!-- of the poisson equation will increase with the number of cycles.
80!-- If the number of cycles is preset by the user, this number will be
81!-- carried out regardless of the accuracy.
82    grid_level_count =   0
83    mgcycles         =   0
84    IF ( mg_cycles == -1 )  THEN
85       maximum_mgcycles = 0
86       residual_norm    = 1.0 
87    ELSE
88       maximum_mgcycles = mg_cycles
89       residual_norm    = 0.0
90    ENDIF
91
92    DO WHILE ( residual_norm > residual_limit  .OR. &
93               mgcycles < maximum_mgcycles )
94
95       CALL next_mg_level( d, p, p3, r)
96       
97!
98!--    Calculate the residual if the user has not preset the number of
99!--    cycles to be performed
100       IF ( maximum_mgcycles == 0 )  THEN
101          CALL resid( d, p, r )
102          maxerror = SUM( r(nzb+1:nzt,nys:nyn,nxl:nxr)**2 )
103#if defined( __parallel )
104          CALL MPI_ALLREDUCE( maxerror, residual_norm, 1, MPI_REAL, MPI_SUM, &
105                              comm2d, ierr)
106#else
107          residual_norm = maxerror
108#endif
109          residual_norm = SQRT( residual_norm )
110       ENDIF
111
112       mgcycles = mgcycles + 1
113
114!
115!--    If the user has not limited the number of cycles, stop the run in case
116!--    of insufficient convergence
117       IF ( mgcycles > 1000  .AND.  mg_cycles == -1 )  THEN
118          IF ( myid == 0 )  THEN
119             PRINT*, '+++ poismg: no sufficient convergence within 1000 cycles'
120          ENDIF
121          CALL local_stop
122       ENDIF
123
124    ENDDO
125
126    DEALLOCATE( p3 )
127
128    CALL cpu_log( log_point_s(29), 'poismg', 'stop' )
129
130 END SUBROUTINE poismg
131
132
133
134 SUBROUTINE resid( f_mg, p_mg, r )
135
136!------------------------------------------------------------------------------!
137! Description:
138! ------------
139! Computes the residual of the perturbation pressure.
140!------------------------------------------------------------------------------!
141
142    USE arrays_3d
143    USE control_parameters
144    USE grid_variables
145    USE indices
146    USE pegrid
147
148    IMPLICIT NONE
149
150    INTEGER ::  i, j, k, l
151
152    REAL, DIMENSION(nzb:nzt_mg(grid_level)+1,                                  &
153                    nys_mg(grid_level)-1:nyn_mg(grid_level)+1,                 &
154                    nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f_mg, p_mg, r
155
156!
157!-- Calculate the residual
158    l = grid_level
159
160!$OMP PARALLEL PRIVATE (i,j,k)
161!$OMP DO
162    DO  i = nxl_mg(l), nxr_mg(l)
163       DO  j = nys_mg(l), nyn_mg(l) 
164          DO  k = nzb+1, nzt_mg(l)
165             r(k,j,i) = f_mg(k,j,i)                                      &
166                        - ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
167                        - ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
168                        - f2_mg(k,l) * p_mg(k+1,j,i)                     &
169                        - f3_mg(k,l) * p_mg(k-1,j,i)                     &
170                        + f1_mg(k,l) * p_mg(k,j,i)
171          ENDDO
172       ENDDO
173    ENDDO
174!$OMP END PARALLEL
175
176!
177!-- Horizontal boundary conditions
178    CALL exchange_horiz( r, 0, 0 )
179
180    IF ( bc_lr /= 'cyclic' )  THEN
181       IF ( inflow_l .OR. outflow_l )  r(:,:,nxl_mg(l)-1) = r(:,:,nxl_mg(l))
182       IF ( inflow_r .OR. outflow_r )  r(:,:,nxr_mg(l)+1) = r(:,:,nxr_mg(l))
183    ENDIF
184
185    IF ( bc_ns /= 'cyclic' )  THEN
186       IF ( inflow_n .OR. outflow_n )  r(:,nyn_mg(l)+1,:) = r(:,nyn_mg(l),:)
187       IF ( inflow_s .OR. outflow_s )  r(:,nys_mg(l)-1,:) = r(:,nys_mg(l),:)
188    ENDIF
189
190!
191!-- Bottom and top boundary conditions
192    IF ( ibc_p_b == 1 )  THEN
193       r(nzb,:,: ) = r(nzb+1,:,:)
194    ELSE
195       r(nzb,:,: ) = 0.0
196    ENDIF
197
198    IF ( ibc_p_t == 1 )  THEN
199       r(nzt_mg(l)+1,:,: ) = r(nzt_mg(l),:,:)
200    ELSE
201       r(nzt_mg(l)+1,:,: ) = 0.0
202    ENDIF
203
204
205 END SUBROUTINE resid
206
207
208
209 SUBROUTINE restrict( f_mg, r )
210
211!------------------------------------------------------------------------------!
212! Description:
213! ------------
214! Interpolates the residual on the next coarser grid with "full weighting"
215! scheme
216!------------------------------------------------------------------------------!
217
218    USE control_parameters
219    USE grid_variables
220    USE indices
221    USE pegrid
222
223    IMPLICIT NONE
224
225    INTEGER ::  i, ic, j, jc, k, kc, l
226
227    REAL, DIMENSION(nzb:nzt_mg(grid_level)+1,                            &
228                    nys_mg(grid_level)-1:nyn_mg(grid_level)+1,           &
229                    nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f_mg
230
231    REAL, DIMENSION(nzb:nzt_mg(grid_level+1)+1,                          &
232                    nys_mg(grid_level+1)-1:nyn_mg(grid_level+1)+1,       &
233                    nxl_mg(grid_level+1)-1:nxr_mg(grid_level+1)+1) ::  r
234
235!
236!-- Interpolate the residual
237    l = grid_level
238
239!$OMP PARALLEL PRIVATE (i,j,k,ic,jc,kc)
240!$OMP DO
241    DO  ic = nxl_mg(l), nxr_mg(l)   
242       i = 2*ic
243       DO  jc = nys_mg(l), nyn_mg(l) 
244          j = 2*jc
245          DO  kc = nzb+1, nzt_mg(l)
246             k = 2*kc-1
247             f_mg(kc,jc,ic) = 1.0 / 64.0 * (                            &
248                              8.0 * r(k,j,i)                            &
249                            + 4.0 * ( r(k,j,i-1)     + r(k,j,i+1)     + &
250                                      r(k,j+1,i)     + r(k,j-1,i)     ) &
251                            + 2.0 * ( r(k,j-1,i-1)   + r(k,j+1,i-1)   + &
252                                      r(k,j-1,i+1)   + r(k,j+1,i+1)   ) &
253                            + 4.0 * r(k-1,j,i)                          &
254                            + 2.0 * ( r(k-1,j,i-1)   + r(k-1,j,i+1)   + &
255                                      r(k-1,j+1,i)   + r(k-1,j-1,i)   ) &
256                            +       ( r(k-1,j-1,i-1) + r(k-1,j+1,i-1) + &
257                                      r(k-1,j-1,i+1) + r(k-1,j+1,i+1) ) &
258                            + 4.0 * r(k+1,j,i)                          &
259                            + 2.0 * ( r(k+1,j,i-1)   + r(k+1,j,i+1)   + &
260                                      r(k+1,j+1,i)   + r(k+1,j-1,i)   ) &
261                            +       ( r(k+1,j-1,i-1) + r(k+1,j+1,i-1) + &
262                                      r(k+1,j-1,i+1) + r(k+1,j+1,i+1) ) &
263                                           ) 
264          ENDDO
265       ENDDO
266    ENDDO
267!$OMP END PARALLEL
268
269!
270!-- Horizontal boundary conditions
271    CALL exchange_horiz( f_mg, 0, 0 )
272
273    IF ( bc_lr /= 'cyclic' )  THEN
274       IF (inflow_l .OR. outflow_l)  f_mg(:,:,nxl_mg(l)-1) = f_mg(:,:,nxl_mg(l))
275       IF (inflow_r .OR. outflow_r)  f_mg(:,:,nxr_mg(l)+1) = f_mg(:,:,nxr_mg(l))
276    ENDIF
277
278    IF ( bc_ns /= 'cyclic' )  THEN
279       IF (inflow_n .OR. outflow_n)  f_mg(:,nyn_mg(l)+1,:) = f_mg(:,nyn_mg(l),:)
280       IF (inflow_s .OR. outflow_s)  f_mg(:,nys_mg(l)-1,:) = f_mg(:,nys_mg(l),:)
281    ENDIF
282
283!
284!-- Bottom and top boundary conditions
285    IF ( ibc_p_b == 1 )  THEN
286       f_mg(nzb,:,: ) = f_mg(nzb+1,:,:)
287    ELSE
288       f_mg(nzb,:,: ) = 0.0
289    ENDIF
290
291    IF ( ibc_p_t == 1 )  THEN
292       f_mg(nzt_mg(l)+1,:,: ) = f_mg(nzt_mg(l),:,:)
293    ELSE
294       f_mg(nzt_mg(l)+1,:,: ) = 0.0
295    ENDIF
296
297
298END SUBROUTINE restrict
299
300
301
302 SUBROUTINE prolong( p, temp )
303
304!------------------------------------------------------------------------------!
305! Description:
306! ------------
307! Interpolates the correction of the perturbation pressure
308! to the next finer grid.
309!------------------------------------------------------------------------------!
310
311    USE control_parameters
312    USE pegrid
313    USE indices
314
315    IMPLICIT NONE
316
317    INTEGER ::  i, j, k, l
318
319    REAL, DIMENSION(nzb:nzt_mg(grid_level-1)+1,                           &
320                    nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1,        &
321                    nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1 ) ::  p
322
323    REAL, DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
324                    nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
325                    nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  temp
326
327
328!
329!-- First, store elements of the coarser grid on the next finer grid
330    l = grid_level
331
332!$OMP PARALLEL PRIVATE (i,j,k)
333!$OMP DO
334    DO  i = nxl_mg(l-1), nxr_mg(l-1)
335       DO  j = nys_mg(l-1), nyn_mg(l-1)
336!CDIR NODEP
337          DO  k = nzb+1, nzt_mg(l-1)
338!
339!--          Points of the coarse grid are directly stored on the next finer
340!--          grid
341             temp(2*k-1,2*j,2*i) = p(k,j,i) 
342!
343!--          Points between two coarse-grid points
344             temp(2*k-1,2*j,2*i+1) = 0.5 * ( p(k,j,i) + p(k,j,i+1) )
345             temp(2*k-1,2*j+1,2*i) = 0.5 * ( p(k,j,i) + p(k,j+1,i) )
346             temp(2*k,2*j,2*i)     = 0.5 * ( p(k,j,i) + p(k+1,j,i) )
347!
348!--          Points in the center of the planes stretched by four points
349!--          of the coarse grid cube
350             temp(2*k-1,2*j+1,2*i+1) = 0.25 * ( p(k,j,i)   + p(k,j,i+1) + &
351                                                p(k,j+1,i) + p(k,j+1,i+1) )
352             temp(2*k,2*j,2*i+1)     = 0.25 * ( p(k,j,i)   + p(k,j,i+1) + &
353                                                p(k+1,j,i) + p(k+1,j,i+1) )
354             temp(2*k,2*j+1,2*i)     = 0.25 * ( p(k,j,i)   + p(k,j+1,i) + &
355                                                p(k+1,j,i) + p(k+1,j+1,i) )
356!
357!--          Points in the middle of coarse grid cube
358             temp(2*k,2*j+1,2*i+1) = 0.125 * ( p(k,j,i)     + p(k,j,i+1)   + &
359                                               p(k,j+1,i)   + p(k,j+1,i+1) + &
360                                               p(k+1,j,i)   + p(k+1,j,i+1) + &
361                                               p(k+1,j+1,i) + p(k+1,j+1,i+1) )
362          ENDDO
363       ENDDO
364    ENDDO
365!$OMP END PARALLEL
366                         
367!
368!-- Horizontal boundary conditions
369    CALL exchange_horiz( temp, 0, 0 )
370
371    IF ( bc_lr /= 'cyclic' )  THEN
372       IF (inflow_l .OR. outflow_l)  temp(:,:,nxl_mg(l)-1) = temp(:,:,nxl_mg(l))
373       IF (inflow_r .OR. outflow_r)  temp(:,:,nxr_mg(l)+1) = temp(:,:,nxr_mg(l))
374    ENDIF
375
376    IF ( bc_ns /= 'cyclic' )  THEN
377       IF (inflow_n .OR. outflow_n)  temp(:,nyn_mg(l)+1,:) = temp(:,nyn_mg(l),:)
378       IF (inflow_s .OR. outflow_s)  temp(:,nys_mg(l)-1,:) = temp(:,nys_mg(l),:)
379    ENDIF
380
381!
382!-- Bottom and top boundary conditions
383    IF ( ibc_p_b == 1 )  THEN
384       temp(nzb,:,: ) = temp(nzb+1,:,:)
385    ELSE
386       temp(nzb,:,: ) = 0.0
387    ENDIF
388
389    IF ( ibc_p_t == 1 )  THEN
390       temp(nzt_mg(l)+1,:,: ) = temp(nzt_mg(l),:,:)
391    ELSE
392       temp(nzt_mg(l)+1,:,: ) = 0.0
393    ENDIF
394
395 
396 END SUBROUTINE prolong
397
398
399 SUBROUTINE redblack( f_mg, p_mg )
400
401!------------------------------------------------------------------------------!
402! Description:
403! ------------
404! Relaxation method for the multigrid scheme. A Gauss-Seidel iteration with
405! 3D-Red-Black decomposition (GS-RB) is used.
406!------------------------------------------------------------------------------!
407
408    USE arrays_3d
409    USE control_parameters
410    USE cpulog
411    USE grid_variables
412    USE indices
413    USE interfaces
414    USE pegrid
415
416    IMPLICIT NONE
417
418    INTEGER :: colour, i, ic, j, jc, jj, k, l, n
419
420    LOGICAL :: unroll
421
422    REAL, DIMENSION(nzb:nzt_mg(grid_level)+1,                                 &
423                    nys_mg(grid_level)-1:nyn_mg(grid_level)+1,                &
424                    nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f_mg, p_mg
425
426
427    l = grid_level
428
429    unroll = ( MOD( nyn_mg(l)-nys_mg(l)+1, 4 ) == 0  .AND. &
430               MOD( nxr_mg(l)-nxl_mg(l)+1, 2 ) == 0 )
431
432    DO  n = 1, ngsrb
433       
434       DO  colour = 1, 2
435
436          IF ( .NOT. unroll )  THEN
437             CALL cpu_log( log_point_s(36), 'redblack_no_unroll', 'start' )
438
439!
440!--          Without unrolling of loops, no cache optimization
441             DO  i = nxl_mg(l), nxr_mg(l), 2
442                DO  j = nys_mg(l) + 2 - colour, nyn_mg(l), 2 
443                   DO  k = nzb+1, nzt_mg(l), 2
444                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
445                                ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
446                              + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
447                              + f2_mg(k,l) * p_mg(k+1,j,i)                     &
448                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
449                                                       )
450                   ENDDO
451                ENDDO
452             ENDDO
453   
454             DO  i = nxl_mg(l)+1, nxr_mg(l), 2
455                DO  j = nys_mg(l) + (colour-1), nyn_mg(l), 2
456                   DO  k = nzb+1, nzt_mg(l), 2 
457                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
458                                ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
459                              + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
460                              + f2_mg(k,l) * p_mg(k+1,j,i)                     &
461                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
462                                                       )
463                   ENDDO
464                ENDDO
465             ENDDO
466 
467             DO  i = nxl_mg(l), nxr_mg(l), 2
468                DO  j = nys_mg(l) + (colour-1), nyn_mg(l), 2
469                   DO  k = nzb+2, nzt_mg(l), 2
470                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
471                                ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
472                              + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
473                              + f2_mg(k,l) * p_mg(k+1,j,i)                     &
474                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
475                                                       )
476                   ENDDO
477                ENDDO
478             ENDDO
479
480             DO  i = nxl_mg(l)+1, nxr_mg(l), 2
481                DO  j = nys_mg(l) + 2 - colour, nyn_mg(l), 2
482                   DO  k = nzb+2, nzt_mg(l), 2
483                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
484                                ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
485                              + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
486                              + f2_mg(k,l) * p_mg(k+1,j,i)                     &
487                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
488                                                       )
489                   ENDDO
490                ENDDO
491             ENDDO
492             CALL cpu_log( log_point_s(36), 'redblack_no_unroll', 'stop' )
493
494          ELSE
495
496!
497!--          Loop unrolling along y, only one i loop for better cache use
498             CALL cpu_log( log_point_s(38), 'redblack_unroll', 'start' )
499             DO  ic = nxl_mg(l), nxr_mg(l), 2
500                DO  jc = nys_mg(l), nyn_mg(l), 4
501                   i  = ic
502                   jj = jc+2-colour
503                   DO  k = nzb+1, nzt_mg(l), 2
504                      j = jj
505                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
506                                ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
507                              + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
508                              + f2_mg(k,l) * p_mg(k+1,j,i)                     &
509                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
510                                                       )
511                      j = jj+2
512                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
513                                ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
514                              + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
515                              + f2_mg(k,l) * p_mg(k+1,j,i)                     &
516                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
517                                                       )
518!                      j = jj+4
519!                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
520!                                ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
521!                              + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
522!                              + f2_mg(k,l) * p_mg(k+1,j,i)                     &
523!                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
524!                                                    )
525!                      j = jj+6
526!                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
527!                                ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
528!                              + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
529!                              + f2_mg(k,l) * p_mg(k+1,j,i)                     &
530!                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
531!                                                       )
532                   ENDDO
533   
534                   i  = ic+1
535                   jj = jc+colour-1
536                   DO  k = nzb+1, nzt_mg(l), 2 
537                      j =jj
538                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
539                                ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
540                              + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
541                              + f2_mg(k,l) * p_mg(k+1,j,i)                     &
542                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
543                                                       )
544                      j = jj+2
545                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
546                                ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
547                              + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
548                              + f2_mg(k,l) * p_mg(k+1,j,i)                     &
549                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
550                                                       )
551!                      j = jj+4
552!                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
553!                                ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
554!                              + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
555!                              + f2_mg(k,l) * p_mg(k+1,j,i)                     &
556!                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
557!                                                       )
558!                      j = jj+6
559!                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
560!                                ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
561!                              + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
562!                              + f2_mg(k,l) * p_mg(k+1,j,i)                     &
563!                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
564!                                                       )
565                   ENDDO
566
567                   i  = ic
568                   jj = jc+colour-1
569                   DO  k = nzb+2, nzt_mg(l), 2
570                      j =jj
571                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
572                                ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
573                              + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
574                              + f2_mg(k,l) * p_mg(k+1,j,i)                     &
575                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
576                                                       )
577                      j = jj+2
578                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
579                                ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
580                              + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
581                              + f2_mg(k,l) * p_mg(k+1,j,i)                     &
582                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
583                                                       )
584!                      j = jj+4
585!                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
586!                                ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
587!                              + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
588!                              + f2_mg(k,l) * p_mg(k+1,j,i)                     &
589!                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
590!                                                       )
591!                      j = jj+6
592!                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
593!                                ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
594!                              + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
595!                              + f2_mg(k,l) * p_mg(k+1,j,i)                     &
596!                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
597!                                                       )
598                   ENDDO
599
600                   i  = ic+1
601                   jj = jc+2-colour
602                   DO  k = nzb+2, nzt_mg(l), 2
603                      j =jj
604                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
605                                ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
606                              + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
607                              + f2_mg(k,l) * p_mg(k+1,j,i)                     &
608                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
609                                                       )
610                      j = jj+2
611                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
612                                ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
613                              + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
614                              + f2_mg(k,l) * p_mg(k+1,j,i)                     &
615                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
616                                                       )
617!                      j = jj+4
618!                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
619!                                ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
620!                              + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
621!                              + f2_mg(k,l) * p_mg(k+1,j,i)                     &
622!                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
623!                                                       )
624!                      j = jj+6
625!                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
626!                                ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
627!                              + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
628!                              + f2_mg(k,l) * p_mg(k+1,j,i)                     &
629!                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
630!                                                       )
631                   ENDDO
632
633                ENDDO
634             ENDDO
635             CALL cpu_log( log_point_s(38), 'redblack_unroll', 'stop' )
636
637          ENDIF
638
639!
640!--       Horizontal boundary conditions
641          CALL exchange_horiz( p_mg, 0, 0 )
642
643          IF ( bc_lr /= 'cyclic' )  THEN
644             IF ( inflow_l .OR. outflow_l )  THEN
645                p_mg(:,:,nxl_mg(l)-1) = p_mg(:,:,nxl_mg(l))
646             ENDIF
647             IF ( inflow_r .OR. outflow_r )  THEN
648                p_mg(:,:,nxr_mg(l)+1) = p_mg(:,:,nxr_mg(l))
649             ENDIF
650          ENDIF
651
652          IF ( bc_ns /= 'cyclic' )  THEN
653             IF ( inflow_n .OR. outflow_n )  THEN
654                p_mg(:,nyn_mg(l)+1,:) = p_mg(:,nyn_mg(l),:)
655             ENDIF
656             IF ( inflow_s .OR. outflow_s )  THEN
657                p_mg(:,nys_mg(l)-1,:) = p_mg(:,nys_mg(l),:)
658             ENDIF
659          ENDIF
660
661!
662!--       Bottom and top boundary conditions
663          IF ( ibc_p_b == 1 )  THEN
664             p_mg(nzb,:,: ) = p_mg(nzb+1,:,:)
665          ELSE
666             p_mg(nzb,:,: ) = 0.0
667          ENDIF
668
669          IF ( ibc_p_t == 1 )  THEN
670             p_mg(nzt_mg(l)+1,:,: ) = p_mg(nzt_mg(l),:,:)
671          ELSE
672             p_mg(nzt_mg(l)+1,:,: ) = 0.0
673          ENDIF
674
675       ENDDO
676
677    ENDDO
678
679
680 END SUBROUTINE redblack
681
682
683
684 SUBROUTINE mg_gather( f2, f2_sub )
685
686    USE control_parameters
687    USE cpulog
688    USE indices
689    USE interfaces
690    USE pegrid
691
692    IMPLICIT NONE
693
694    INTEGER ::  n, nwords, sender
695
696    REAL, DIMENSION(nzb:nzt_mg(grid_level)+1,                            &
697                    nys_mg(grid_level)-1:nyn_mg(grid_level)+1,           &
698                    nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f2
699
700    REAL, DIMENSION(nzb:mg_loc_ind(5,myid)+1,                            &
701                    mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1,           &
702                    mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) ::  f2_sub
703
704!
705!-- Find out the number of array elements of the subdomain array
706    nwords = SIZE( f2_sub )
707
708#if defined( __parallel )
709    CALL cpu_log( log_point_s(34), 'mg_gather', 'start' )
710
711    IF ( myid == 0 )  THEN
712!
713!--    Store the local subdomain array on the total array
714       f2(:,mg_loc_ind(3,0)-1:mg_loc_ind(4,0)+1, &
715            mg_loc_ind(1,0)-1:mg_loc_ind(2,0)+1) = f2_sub
716
717!
718!--    Receive the subdomain arrays from all other PEs and store them on the
719!--    total array
720       DO  n = 1, numprocs-1
721!
722!--       Receive the arrays in arbitrary order from the PEs.
723          CALL MPI_RECV( f2_sub(nzb,mg_loc_ind(3,0)-1,mg_loc_ind(1,0)-1),     &
724                         nwords, MPI_REAL, MPI_ANY_SOURCE, 1, comm2d, status, &
725                         ierr )
726          sender = status(MPI_SOURCE)
727          f2(:,mg_loc_ind(3,sender)-1:mg_loc_ind(4,sender)+1, &
728               mg_loc_ind(1,sender)-1:mg_loc_ind(2,sender)+1) = f2_sub
729       ENDDO
730
731    ELSE
732!
733!--    Send subdomain array to PE0
734       CALL MPI_SEND( f2_sub(nzb,mg_loc_ind(3,myid)-1,mg_loc_ind(1,myid)-1), &
735                      nwords, MPI_REAL, 0, 1, comm2d, ierr )
736    ENDIF
737
738    CALL cpu_log( log_point_s(34), 'mg_gather', 'stop' )
739#endif
740   
741 END SUBROUTINE mg_gather
742
743
744
745 SUBROUTINE mg_scatter( p2, p2_sub )
746!
747!-- TODO: It may be possible to improve the speed of this routine by using
748!-- non-blocking communication
749
750    USE control_parameters
751    USE cpulog
752    USE indices
753    USE interfaces
754    USE pegrid
755
756    IMPLICIT NONE
757
758    INTEGER ::  n, nwords, sender
759
760    REAL, DIMENSION(nzb:nzt_mg(grid_level-1)+1,                            &
761                    nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1,         &
762                    nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) ::  p2
763
764    REAL, DIMENSION(nzb:mg_loc_ind(5,myid)+1,                              &
765                    mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1,             &
766                    mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) ::  p2_sub
767
768!
769!-- Find out the number of array elements of the subdomain array
770    nwords = SIZE( p2_sub )
771
772#if defined( __parallel )
773    CALL cpu_log( log_point_s(35), 'mg_scatter', 'start' )
774
775    IF ( myid == 0 )  THEN
776!
777!--    Scatter the subdomain arrays to the other PEs by blocking
778!--    communication
779       DO  n = 1, numprocs-1
780
781          p2_sub = p2(:,mg_loc_ind(3,n)-1:mg_loc_ind(4,n)+1, &
782                        mg_loc_ind(1,n)-1:mg_loc_ind(2,n)+1)
783
784          CALL MPI_SEND( p2_sub(nzb,mg_loc_ind(3,0)-1,mg_loc_ind(1,0)-1), &
785                         nwords, MPI_REAL, n, 1, comm2d, ierr )
786
787       ENDDO
788
789!
790!--    Store data from the total array to the local subdomain array
791       p2_sub = p2(:,mg_loc_ind(3,0)-1:mg_loc_ind(4,0)+1, &
792                     mg_loc_ind(1,0)-1:mg_loc_ind(2,0)+1)
793
794    ELSE
795!
796!--    Receive subdomain array from PE0
797       CALL MPI_RECV( p2_sub(nzb,mg_loc_ind(3,myid)-1,mg_loc_ind(1,myid)-1), &
798                      nwords, MPI_REAL, 0, 1, comm2d, status, ierr )
799
800    ENDIF
801
802    CALL cpu_log( log_point_s(35), 'mg_scatter', 'stop' )
803#endif
804   
805 END SUBROUTINE mg_scatter
806
807
808
809 RECURSIVE SUBROUTINE next_mg_level( f_mg, p_mg, p3, r )
810
811!------------------------------------------------------------------------------!
812! Description:
813! ------------
814! This is where the multigrid technique takes place. V- and W- Cycle are
815! implemented and steered by the parameter "gamma". Parameter "nue" determines
816! the convergence of the multigrid iterative solution. There are nue times
817! RB-GS iterations. It should be set to "1" or "2", considering the time effort
818! one would like to invest. Last choice shows a very good converging factor,
819! but leads to an increase in computing time.
820!------------------------------------------------------------------------------!
821
822    USE arrays_3d
823    USE control_parameters
824    USE grid_variables
825    USE indices
826    USE pegrid
827
828    IMPLICIT NONE
829
830    INTEGER ::  i, j, k, nxl_mg_save, nxr_mg_save, nyn_mg_save, nys_mg_save, &
831                nzt_mg_save
832
833    LOGICAL ::  restore_boundary_lr_on_pe0, restore_boundary_ns_on_pe0
834
835    REAL, DIMENSION(nzb:nzt_mg(grid_level)+1,                                  &
836                 nys_mg(grid_level)-1:nyn_mg(grid_level)+1,                    &
837                 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg, p_mg, p3, r
838
839    REAL, DIMENSION(:,:,:), ALLOCATABLE ::  f2, f2_sub, p2, p2_sub
840
841!
842!-- Restriction to the coarsest grid
843 10 IF ( grid_level == 1 )  THEN
844
845!
846!--    Solution on the coarsest grid. Double the number of Gauss-Seidel
847!--    iterations in order to get a more accurate solution.
848       ngsrb = 2 * ngsrb
849       CALL redblack( f_mg, p_mg )
850       ngsrb = ngsrb / 2
851
852    ELSEIF ( grid_level /= 1 )  THEN
853
854       grid_level_count(grid_level) = grid_level_count(grid_level) + 1
855
856!
857!--    Solution on the actual grid level
858       CALL redblack( f_mg, p_mg )
859
860!
861!--    Determination of the actual residual
862       CALL resid( f_mg, p_mg, r )
863
864!
865!--    Restriction of the residual (finer grid values!) to the next coarser
866!--    grid. Therefore, the grid level has to be decremented now. nxl..nzt have
867!--    to be set to the coarse grid values, because these variables are needed
868!--    for the exchange of ghost points in routine exchange_horiz
869       grid_level = grid_level - 1
870       nxl = nxl_mg(grid_level)
871       nxr = nxr_mg(grid_level)
872       nys = nys_mg(grid_level)
873       nyn = nyn_mg(grid_level)
874       nzt = nzt_mg(grid_level)
875
876       ALLOCATE( f2(nzb:nzt_mg(grid_level)+1,                    &
877                    nys_mg(grid_level)-1:nyn_mg(grid_level)+1,   &
878                    nxl_mg(grid_level)-1:nxr_mg(grid_level)+1),  &
879                 p2(nzb:nzt_mg(grid_level)+1,                    &
880                    nys_mg(grid_level)-1:nyn_mg(grid_level)+1,   &
881                    nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) )
882
883       IF ( grid_level == mg_switch_to_pe0_level )  THEN
884!          print*, 'myid=',myid, ' restrict and switch to PE0. level=', grid_level
885!
886!--       From this level on, calculations are done on PE0 only.
887!--       First, carry out restriction on the subdomain.
888!--       Therefore, indices of the level have to be changed to subdomain values
889!--       in between (otherwise, the restrict routine would expect
890!--       the gathered array)
891          nxl_mg_save = nxl_mg(grid_level)
892          nxr_mg_save = nxr_mg(grid_level)
893          nys_mg_save = nys_mg(grid_level)
894          nyn_mg_save = nyn_mg(grid_level)
895          nzt_mg_save = nzt_mg(grid_level)
896          nxl_mg(grid_level) = mg_loc_ind(1,myid)
897          nxr_mg(grid_level) = mg_loc_ind(2,myid)
898          nys_mg(grid_level) = mg_loc_ind(3,myid)
899          nyn_mg(grid_level) = mg_loc_ind(4,myid)
900          nzt_mg(grid_level) = mg_loc_ind(5,myid)
901          nxl = mg_loc_ind(1,myid)
902          nxr = mg_loc_ind(2,myid)
903          nys = mg_loc_ind(3,myid)
904          nyn = mg_loc_ind(4,myid)
905          nzt = mg_loc_ind(5,myid)
906
907          ALLOCATE( f2_sub(nzb:nzt_mg(grid_level)+1,                    &
908                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,   &
909                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) )
910
911          CALL restrict( f2_sub, r )
912
913!
914!--       Restore the correct indices of this level
915          nxl_mg(grid_level) = nxl_mg_save
916          nxr_mg(grid_level) = nxr_mg_save
917          nys_mg(grid_level) = nys_mg_save
918          nyn_mg(grid_level) = nyn_mg_save
919          nzt_mg(grid_level) = nzt_mg_save
920          nxl = nxl_mg(grid_level)
921          nxr = nxr_mg(grid_level)
922          nys = nys_mg(grid_level)
923          nyn = nyn_mg(grid_level)
924          nzt = nzt_mg(grid_level)
925
926!
927!--       Gather all arrays from the subdomains on PE0
928          CALL mg_gather( f2, f2_sub )
929
930!
931!--       Set switch for routine exchange_horiz, that no ghostpoint exchange
932!--       has to be carried out from now on
933          mg_switch_to_pe0 = .TRUE.
934
935!
936!--       In case of non-cyclic lateral boundary conditions, both in- and
937!--       outflow conditions have to be used on PE0 after the switch, because
938!--       it then contains the total domain. Due to the virtual processor
939!--       grid, before the switch, PE0 can have in-/outflow at the left
940!--       and south wall only (or on opposite walls in case of a 1d
941!--       decomposition).
942          restore_boundary_lr_on_pe0 = .FALSE.
943          restore_boundary_ns_on_pe0 = .FALSE.
944          IF ( myid == 0 )  THEN
945             IF ( inflow_l  .AND.  .NOT. outflow_r )  THEN
946                outflow_r = .TRUE.
947                restore_boundary_lr_on_pe0 = .TRUE.
948             ENDIF
949             IF ( outflow_l  .AND.  .NOT. inflow_r )  THEN
950                inflow_r  = .TRUE.
951                restore_boundary_lr_on_pe0 = .TRUE.
952             ENDIF
953             IF ( inflow_s  .AND.  .NOT. outflow_n )  THEN
954                outflow_n = .TRUE.
955                restore_boundary_ns_on_pe0 = .TRUE.
956             ENDIF
957             IF ( outflow_s  .AND.  .NOT. inflow_n )  THEN
958                inflow_n  = .TRUE.
959                restore_boundary_ns_on_pe0 = .TRUE.
960             ENDIF
961          ENDIF
962
963          DEALLOCATE( f2_sub )
964
965       ELSE
966
967          CALL restrict( f2, r )
968
969       ENDIF
970       p2 = 0.0
971
972!
973!--    Repeat the same procedure till the coarsest grid is reached
974       IF ( myid == 0  .OR.  grid_level > mg_switch_to_pe0_level )  THEN
975          CALL next_mg_level( f2, p2, p3, r )
976       ENDIF
977
978    ENDIF
979
980!
981!-- Now follows the prolongation
982    IF ( grid_level >= 2 )  THEN
983
984!
985!--    Grid level has to be incremented on the PEs where next_mg_level
986!--    has not been called before (normally it is incremented at the end
987!--    of next_mg_level)
988       IF ( myid /= 0  .AND.  grid_level == mg_switch_to_pe0_level )  THEN
989          grid_level = grid_level + 1
990          nxl = nxl_mg(grid_level)
991          nxr = nxr_mg(grid_level)
992          nys = nys_mg(grid_level)
993          nyn = nyn_mg(grid_level)
994          nzt = nzt_mg(grid_level)
995       ENDIF
996
997!
998!--    Prolongation of the new residual. The values are transferred
999!--    from the coarse to the next finer grid.
1000       IF ( grid_level == mg_switch_to_pe0_level+1 )  THEN
1001!
1002!--       At this level, the new residual first has to be scattered from
1003!--       PE0 to the other PEs
1004          ALLOCATE( p2_sub(nzb:mg_loc_ind(5,myid)+1,             &
1005                    mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1,   &
1006                    mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) )
1007
1008          CALL mg_scatter( p2, p2_sub )
1009
1010!
1011!--       Therefore, indices of the previous level have to be changed to
1012!--       subdomain values in between (otherwise, the prolong routine would
1013!--       expect the gathered array)
1014          nxl_mg_save = nxl_mg(grid_level-1)
1015          nxr_mg_save = nxr_mg(grid_level-1)
1016          nys_mg_save = nys_mg(grid_level-1)
1017          nyn_mg_save = nyn_mg(grid_level-1)
1018          nzt_mg_save = nzt_mg(grid_level-1)
1019          nxl_mg(grid_level-1) = mg_loc_ind(1,myid)
1020          nxr_mg(grid_level-1) = mg_loc_ind(2,myid)
1021          nys_mg(grid_level-1) = mg_loc_ind(3,myid)
1022          nyn_mg(grid_level-1) = mg_loc_ind(4,myid)
1023          nzt_mg(grid_level-1) = mg_loc_ind(5,myid)
1024
1025!
1026!--       Set switch for routine exchange_horiz, that ghostpoint exchange
1027!--       has to be carried again out from now on
1028          mg_switch_to_pe0 = .FALSE.
1029
1030!
1031!--       In case of non-cyclic lateral boundary conditions, restore the
1032!--       in-/outflow conditions on PE0
1033          IF ( myid == 0 )  THEN
1034             IF ( restore_boundary_lr_on_pe0 )  THEN
1035                IF ( inflow_l  )  outflow_r = .FALSE.
1036                IF ( outflow_l )  inflow_r  = .FALSE.
1037             ENDIF
1038             IF ( restore_boundary_ns_on_pe0 )  THEN
1039                IF ( inflow_s  )  outflow_n = .FALSE.
1040                IF ( outflow_s )  inflow_n  = .FALSE.
1041             ENDIF
1042          ENDIF
1043
1044          CALL prolong( p2_sub, p3 )
1045
1046!
1047!--       Restore the correct indices of the previous level
1048          nxl_mg(grid_level-1) = nxl_mg_save
1049          nxr_mg(grid_level-1) = nxr_mg_save
1050          nys_mg(grid_level-1) = nys_mg_save
1051          nyn_mg(grid_level-1) = nyn_mg_save
1052          nzt_mg(grid_level-1) = nzt_mg_save
1053
1054          DEALLOCATE( p2_sub )
1055
1056       ELSE
1057
1058          CALL prolong( p2, p3 )
1059
1060       ENDIF
1061
1062!
1063!--    Temporary arrays for the actual grid are not needed any more
1064       DEALLOCATE( p2, f2 )
1065
1066!
1067!--    Computation of the new pressure correction. Therefore,
1068!--    values from prior grids are added up automatically stage by stage.
1069       DO  i = nxl_mg(grid_level)-1, nxr_mg(grid_level)+1
1070          DO  j = nys_mg(grid_level)-1, nyn_mg(grid_level)+1
1071             DO  k = nzb, nzt_mg(grid_level)+1 
1072                p_mg(k,j,i) = p_mg(k,j,i) + p3(k,j,i)
1073             ENDDO
1074          ENDDO
1075       ENDDO
1076
1077!
1078!--    Relaxation of the new solution
1079       CALL redblack( f_mg, p_mg )
1080
1081    ENDIF
1082
1083!
1084!-- The following few lines serve the steering of the multigrid scheme
1085    IF ( grid_level == maximum_grid_level )  THEN
1086
1087       GOTO 20
1088
1089    ELSEIF ( grid_level /= maximum_grid_level  .AND.  grid_level /= 1  .AND. &
1090             grid_level_count(grid_level) /= gamma_mg )  THEN
1091
1092       GOTO 10
1093
1094    ENDIF
1095
1096!
1097!-- Reset counter for the next call of poismg
1098    grid_level_count(grid_level) = 0
1099
1100!
1101!-- Continue with the next finer level. nxl..nzt have to be
1102!-- set to the finer grid values, because these variables are needed for the
1103!-- exchange of ghost points in routine exchange_horiz
1104    grid_level = grid_level + 1
1105    nxl = nxl_mg(grid_level)
1106    nxr = nxr_mg(grid_level)
1107    nys = nys_mg(grid_level)
1108    nyn = nyn_mg(grid_level)
1109    nzt = nzt_mg(grid_level)
1110
1111 20 CONTINUE
1112
1113 END SUBROUTINE next_mg_level
Note: See TracBrowser for help on using the repository browser.