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

Last change on this file since 77 was 77, checked in by raasch, 14 years ago

New:
---

particle reflection from vertical walls implemented, particle SGS model adjusted to walls

Wall functions for vertical walls now include diabatic conditions. New subroutines wall_fluxes, wall_fluxes_e. New 4D-array rif_wall.

new d3par-parameter netcdf_64bit_3d to switch on 64bit offset only for 3D files

new d3par-parameter dt_max to define the maximum value for the allowed timestep

new inipar-parameter loop_optimization to control the loop optimization method

new inipar-parameter pt_refrence. If given, this value is used as the reference that in buoyancy terms (otherwise, the instantaneous horizontally averaged temperature is used).

new user interface user_advec_particles

new initializing action "by_user" calls user_init_3d_model and allows the initial setting of all 3d arrays

topography height informations are stored on arrays zu_s_inner and zw_w_inner and output to the 2d/3d NetCDF files

samples added to the user interface which show how to add user-define time series quantities.

calculation/output of precipitation amount, precipitation rate and z0 (by setting "pra*", "prr*", "z0*" with data_output). The time interval on which the precipitation amount is defined is set by new d3par-parameter precipitation_amount_interval

unit 9 opened for debug output (file DEBUG_<pe#>)

Makefile, advec_particles, average_3d_data, buoyancy, calc_precipitation, check_open, check_parameters, data_output_2d, diffusion_e, diffusion_u, diffusion_v, diffusion_w, diffusivities, header, impact_of_latent_heat, init_particles, init_3d_model, modules, netcdf, parin, production_e, read_var_list, read_3d_binary, sum_up_3d_data, user_interface, write_var_list, write_3d_binary

New: wall_fluxes

Changed:


General revision of non-cyclic horizontal boundary conditions: radiation boundary conditions are now used instead of Neumann conditions at the outflow (calculation needs velocity values for t-dt, which are stored on new arrays u_m_l, u_m_r, etc.), calculation of mean outflow is not needed any more, volume flow control is added for the outflow boundary (currently only for the north boundary!!), additional gridpoints along x and y (uxrp, vynp) are not needed any more, routine "boundary_conds" now operates on timelevel t+dt and is not split in two parts (main, uvw_outflow) any more, Neumann boundary conditions at inflow/outflow in case of non-cyclic boundary conditions for all 2d-arrays that are handled by exchange_horiz_2d

The FFT-method for solving the Poisson-equation is now working with Neumann boundary conditions both at the bottom and the top. This requires adjustments of the tridiagonal coefficients and subtracting the horizontally averaged mean from the vertical velocity field.

+age_m in particle_type

Particles-package is now part of the default code ("-p particles" is not needed any more).

Move call of user_actions( 'after_integration' ) below increment of times
and counters. user_actions is now called for each statistic region and has as an argument the number of the respective region (sr)

d3par-parameter data_output_ts removed. Timeseries output for "profil" removed. Timeseries are now switched on by dt_dots. Timeseries data is collected in flow_statistics.

Initial velocities at nzb+1 are regarded for volume flow control in case they have been set zero before (to avoid small timesteps); see new internal parameters u/v_nzb_p1_for_vfc.

q is not allowed to become negative (prognostic_equations).

poisfft_init is only called if fft-solver is switched on (init_pegrid).

d3par-parameter moisture renamed to humidity.

Subversion global revision number is read from mrun and added to the run description header and to the run control (_rc) file.

vtk directives removed from main program.

The uitility routine interpret_config reads PALM environment variables from NAMELIST instead using the system call GETENV.

advec_u_pw, advec_u_up, advec_v_pw, advec_v_up, asselin_filter, check_parameters, coriolis, data_output_dvrp, data_output_ptseries, data_output_ts, data_output_2d, data_output_3d, diffusion_u, diffusion_v, exchange_horiz, exchange_horiz_2d, flow_statistics, header, init_grid, init_particles, init_pegrid, init_rankine, init_pt_anomaly, init_1d_model, init_3d_model, modules, palm, package_parin, parin, poisfft, poismg, prandtl_fluxes, pres, production_e, prognostic_equations, read_var_list, read_3d_binary, sor, swap_timelevel, time_integration, write_var_list, write_3d_binary

Errors:


Bugfix: preset of tendencies te_em, te_um, te_vm in init_1d_model

Bugfix in sample for reading user defined data from restart file (user_init)

Bugfix in setting diffusivities for cases with the outflow damping layer extending over more than one subdomain (init_3d_model)

Check for possible negative humidities in the initial humidity profile.

in Makefile, default suffixes removed from the suffix list to avoid calling of m2c in
# case of .mod files

Makefile
check_parameters, init_1d_model, init_3d_model, user_interface

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