source: palm/trunk/SOURCE/poismg_mod.f90 @ 4648

Last change on this file since 4648 was 4457, checked in by raasch, 5 years ago

ghost point exchange modularized, bugfix for wrong 2d-exchange

  • Property svn:keywords set to Id
File size: 85.9 KB
RevLine 
[1931]1!> @file poismg.f90
[2000]2!------------------------------------------------------------------------------!
[2696]3! This file is part of the PALM model system.
[1575]4!
[2000]5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
[1575]9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
[4360]17! Copyright 1997-2020 Leibniz Universitaet Hannover
[2000]18!------------------------------------------------------------------------------!
[1575]19!
20! Current revisions:
21! -----------------
[1610]22!
[3183]23!
[1576]24! Former revisions:
25! -----------------
26! $Id: poismg_mod.f90 4457 2020-03-11 14:20:43Z raasch $
[4457]27! use statement for exchange horiz added
28!
29! 4432 2020-02-28 07:43:21Z raasch
[4432]30! bugfix for previous revision (vector directive was changed by mistake)
31!
32! 4429 2020-02-27 15:24:30Z raasch
[4429]33! statement added to avoid compile error due to unused dummy argument
34! bugfix: cpp-directives added for serial mode
35!
36! 4360 2020-01-07 11:25:50Z suehring
[4182]37! Corrected "Former revisions" section
38!
39! 3725 2019-02-07 10:11:02Z raasch
[3725]40! unused subroutine removed
41!
42! 3655 2019-01-07 16:51:22Z knoop
[3534]43! unnecessary check eliminated
[3182]44!
[4182]45! Following optimisations have been made:
46! - vectorisation (for Intel-CPUs) of the red-black algorithm by resorting
47!   array elements with even and odd indices
48! - explicit boundary conditions for building walls removed (solver is
49!   running through the buildings
50! - reduced data transfer in case of ghost point exchange, because only
51!   "red" or "black" data points need to be exchanged. This is not applied
52!   for coarser grid levels, since for then the transfer time is latency bound
53!
54!
[1575]55! Description:
56! ------------
[1682]57!> Solves the Poisson equation for the perturbation pressure with a multigrid
58!> V- or W-Cycle scheme.
59!>
60!> This multigrid method was originally developed for PALM by Joerg Uhlenbrock,
61!> September 2000 - July 2001. It has been optimised for speed by Klaus
62!> Ketelsen in November 2014.
63!>
64!> @attention Loop unrolling and cache optimization in SOR-Red/Black method
65!>            still does not give the expected speedup!
66!>
67!> @todo Further work required.
[1575]68!------------------------------------------------------------------------------!
[1682]69 MODULE poismg_mod
70 
[2696]71    USE control_parameters,                                                    &
[3182]72        ONLY:  bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r,                 &
73               bc_dirichlet_s, bc_radiation_l, bc_radiation_n, bc_radiation_r, &
74               bc_radiation_s, grid_level, nesting_offline
[1575]75
76    USE cpulog,                                                                &
77        ONLY:  cpu_log, log_point_s
[1609]78
[4457]79    USE exchange_horiz_mod,                                                    &
80        ONLY:  exchange_horiz
81
[1575]82    USE kinds
[1609]83
[1575]84    USE pegrid
85
86    PRIVATE
87
[1682]88    INTEGER, SAVE                             ::  ind_even_odd    !< border index between even and odd k index
89    INTEGER, DIMENSION(:), SAVE, ALLOCATABLE  ::  even_odd_level  !< stores ind_even_odd for all MG levels
[1575]90
[1682]91    REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE ::  f1_mg_b, f2_mg_b, f3_mg_b  !< blocked version of f1_mg ...
[1575]92
[1931]93    INTERFACE poismg
94       MODULE PROCEDURE poismg
95    END INTERFACE poismg
[1575]96
97    INTERFACE sort_k_to_even_odd_blocks
98       MODULE PROCEDURE sort_k_to_even_odd_blocks
[3725]99!       MODULE PROCEDURE sort_k_to_even_odd_blocks_int
[1575]100       MODULE PROCEDURE sort_k_to_even_odd_blocks_1d
101    END INTERFACE sort_k_to_even_odd_blocks
102
[1931]103    PUBLIC poismg
[1575]104
105 CONTAINS
106
[1682]107!------------------------------------------------------------------------------!
108! Description:
109! ------------
110!> Solves the Poisson equation for the perturbation pressure with a multigrid
111!> V- or W-Cycle scheme.
112!------------------------------------------------------------------------------!
[1931]113    SUBROUTINE poismg( r )
[1575]114
115       USE arrays_3d,                                                          &
116           ONLY:  d, p_loc
117
118       USE control_parameters,                                                 &
[2939]119           ONLY:  bc_lr_cyc, bc_ns_cyc, gathered_size, grid_level,             &
120                  grid_level_count, ibc_p_t,                                   &
[1575]121                  maximum_grid_level, message_string, mgcycles, mg_cycles,     &
122                  mg_switch_to_pe0_level, residual_limit, subdomain_size
123
124       USE cpulog,                                                             &
125           ONLY:  cpu_log, log_point_s
126
127       USE indices,                                                            &
128           ONLY:  nxl, nxlg, nxl_mg, nxr, nxrg, nxr_mg, nys, nysg, nys_mg, nyn,&
129                  nyng, nyn_mg, nzb, nzt, nzt_mg
130
131       IMPLICIT NONE
132
[1682]133       REAL(wp) ::  maxerror          !<
134       REAL(wp) ::  maximum_mgcycles  !<
135       REAL(wp) ::  residual_norm     !<
[1575]136
[1682]137       REAL(wp), DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) ::  r  !<
[1575]138
[1682]139       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  p3  !<
[1575]140
141
[1931]142       CALL cpu_log( log_point_s(29), 'poismg', 'start' )
[1575]143!
144!--    Initialize arrays and variables used in this subroutine
145
146!--    If the number of grid points of the gathered grid, which is collected
147!--    on PE0, is larger than the number of grid points of an PE, than array
148!--    p3 will be enlarged.
149       IF ( gathered_size > subdomain_size )  THEN
150          ALLOCATE( p3(nzb:nzt_mg(mg_switch_to_pe0_level)+1,nys_mg(            &
151                    mg_switch_to_pe0_level)-1:nyn_mg(mg_switch_to_pe0_level)+1,&
152                    nxl_mg(mg_switch_to_pe0_level)-1:nxr_mg(                   &
153                    mg_switch_to_pe0_level)+1) )
154       ELSE
155          ALLOCATE ( p3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
156       ENDIF
157
158       p3 = 0.0_wp
159
160 
161!
162!--    Ghost boundaries have to be added to divergence array.
163!--    Exchange routine needs to know the grid level!
164       grid_level = maximum_grid_level
165       CALL exchange_horiz( d, 1)
[2939]166!
167!--    Set bottom and top boundary conditions
[1575]168       d(nzb,:,:) = d(nzb+1,:,:)
[2939]169       IF ( ibc_p_t == 1 )  d(nzt+1,:,: ) = d(nzt,:,:)
[1575]170!
[2939]171!--    Set lateral boundary conditions in non-cyclic case
172       IF ( .NOT. bc_lr_cyc )  THEN
[3182]173          IF ( bc_dirichlet_l  .OR.  bc_radiation_l )                          &
[2939]174             d(:,:,nxl-1) = d(:,:,nxl)
[3182]175          IF ( bc_dirichlet_r  .OR.  bc_radiation_r )                          &
[2939]176             d(:,:,nxr+1) = d(:,:,nxr)
177       ENDIF
178       IF ( .NOT. bc_ns_cyc )  THEN
[3182]179          IF ( bc_dirichlet_n  .OR.  bc_radiation_n )                          &
[2939]180             d(:,nyn+1,:) = d(:,nyn,:)
[3182]181          IF ( bc_dirichlet_s  .OR.  bc_radiation_s )                          &
[2939]182             d(:,nys-1,:) = d(:,nys,:)
183       ENDIF
184!
[1575]185!--    Initiation of the multigrid scheme. Does n cycles until the
186!--    residual is smaller than the given limit. The accuracy of the solution
187!--    of the poisson equation will increase with the number of cycles.
188!--    If the number of cycles is preset by the user, this number will be
189!--    carried out regardless of the accuracy.
190       grid_level_count =  0
191       mgcycles         =  0
192       IF ( mg_cycles == -1 )  THEN
193          maximum_mgcycles = 0
194          residual_norm    = 1.0_wp
195       ELSE
196          maximum_mgcycles = mg_cycles
197          residual_norm    = 0.0_wp
198       ENDIF
199
200!
201!--    Initial settings for sorting k-dimension from sequential order (alternate
202!--    even/odd) into blocks of even and odd or vice versa
203       CALL init_even_odd_blocks
204
205!
206!--    Sort input arrays in even/odd blocks along k-dimension
207       CALL sort_k_to_even_odd_blocks( d, grid_level )
208       CALL sort_k_to_even_odd_blocks( p_loc, grid_level )
209
210!
211!--    The complete multigrid cycles are running in block mode, i.e. over
212!--    seperate data blocks of even and odd indices
213       DO WHILE ( residual_norm > residual_limit  .OR. &
214                  mgcycles < maximum_mgcycles )
215 
[1931]216          CALL next_mg_level( d, p_loc, p3, r)
[1575]217
218!
219!--       Calculate the residual if the user has not preset the number of
220!--       cycles to be performed
221          IF ( maximum_mgcycles == 0 )  THEN
[1931]222             CALL resid( d, p_loc, r )
[1575]223             maxerror = SUM( r(nzb+1:nzt,nys:nyn,nxl:nxr)**2 )
224
225#if defined( __parallel )
226             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
227                CALL MPI_ALLREDUCE( maxerror, residual_norm, 1, MPI_REAL,      &
228                                    MPI_SUM, comm2d, ierr)
229#else
230                residual_norm = maxerror
231#endif
232             residual_norm = SQRT( residual_norm )
233          ENDIF
234
235          mgcycles = mgcycles + 1
236
237!
238!--       If the user has not limited the number of cycles, stop the run in case
239!--       of insufficient convergence
240          IF ( mgcycles > 1000  .AND.  mg_cycles == -1 )  THEN
241             message_string = 'no sufficient convergence within 1000 cycles'
[1931]242             CALL message( 'poismg', 'PA0283', 1, 2, 0, 6, 0 )
[1575]243          ENDIF
244
245       ENDDO
246
247       DEALLOCATE( p3 )
248!
249!--    Result has to be sorted back from even/odd blocks to sequential order
250       CALL sort_k_to_sequential( p_loc )
251!
252!--    Unset the grid level. Variable is used to determine the MPI datatypes for
253!--    ghost point exchange
254       grid_level = 0
255
[1931]256       CALL cpu_log( log_point_s(29), 'poismg', 'stop' )
[1575]257
[1931]258    END SUBROUTINE poismg
[1575]259
260
261!------------------------------------------------------------------------------!
262! Description:
263! ------------
[1682]264!> Computes the residual of the perturbation pressure.
[1575]265!------------------------------------------------------------------------------!
[1931]266    SUBROUTINE resid( f_mg, p_mg, r )
[1575]267
[1682]268
[1575]269       USE arrays_3d,                                                          &
[3241]270           ONLY:  rho_air_mg
[1575]271
272       USE control_parameters,                                                 &
[2696]273           ONLY:  bc_lr_cyc, bc_ns_cyc, ibc_p_b, ibc_p_t
[1575]274       USE grid_variables,                                                     &
275           ONLY:  ddx2_mg, ddy2_mg
276
277       USE indices,                                                            &
[1898]278           ONLY:  nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
[1575]279
280       IMPLICIT NONE
281
[1898]282       INTEGER(iwp) ::  i        !< index variable along x
283       INTEGER(iwp) ::  j        !< index variable along y
284       INTEGER(iwp) ::  k        !< index variable along z
285       INTEGER(iwp) ::  l        !< index indicating grid level
286       INTEGER(iwp) ::  km1      !< index variable along z dimension (k-1)
287       INTEGER(iwp) ::  kp1      !< index variable along z dimension (k+1)
[1575]288
289       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
290                       nys_mg(grid_level)-1:nyn_mg(grid_level)+1,              &
[1898]291                       nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f_mg  !< velocity divergence
[1575]292       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
293                       nys_mg(grid_level)-1:nyn_mg(grid_level)+1,              &
[1898]294                       nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  p_mg  !< perturbation pressure
[1575]295       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
296                       nys_mg(grid_level)-1:nyn_mg(grid_level)+1,              &
[1898]297                       nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  r     !< residuum of perturbation pressure
[1575]298
299!
300!--    Calculate the residual
301       l = grid_level
302
303       CALL cpu_log( log_point_s(53), 'resid', 'start' )
[1898]304       !$OMP PARALLEL PRIVATE (i,j,k,km1,kp1)
305       !$OMP DO
306       DO  i = nxl_mg(l), nxr_mg(l)
307          DO  j = nys_mg(l), nyn_mg(l)
[1575]308                !DIR$ IVDEP
[1898]309             DO k = ind_even_odd+1, nzt_mg(l)
310                km1 = k-ind_even_odd-1
311                kp1 = k-ind_even_odd
312                r(k,j,i) = f_mg(k,j,i)                                         &
[2232]313                      - rho_air_mg(k,l) * ddx2_mg(l) *                         &
[1898]314                      ( p_mg(k,j,i+1) +  p_mg(k,j,i-1)  )                      &
[2232]315                      - rho_air_mg(k,l) * ddy2_mg(l) *                         &
[1898]316                      ( p_mg(k,j+1,i) + p_mg(k,j-1,i)  )                       &
317                      - f2_mg_b(k,l) * p_mg(kp1,j,i)                           &
318                      - f3_mg_b(k,l) * p_mg(km1,j,i)                           &
[1575]319                      + f1_mg_b(k,l) * p_mg(k,j,i)
[1898]320             ENDDO
321             !DIR$ IVDEP
322             DO k = nzb+1, ind_even_odd
323                km1 = k+ind_even_odd
324                kp1 = k+ind_even_odd+1
325                r(k,j,i) = f_mg(k,j,i)                                         &
[2232]326                      - rho_air_mg(k,l) * ddx2_mg(l) *                         &
[1898]327                      ( p_mg(k,j,i+1) +  p_mg(k,j,i-1)  )                      &
[2232]328                      - rho_air_mg(k,l) * ddy2_mg(l) *                         &
[1898]329                      ( p_mg(k,j+1,i) + p_mg(k,j-1,i)  )                       &
330                      - f2_mg_b(k,l) * p_mg(kp1,j,i)                           &
331                      - f3_mg_b(k,l) * p_mg(km1,j,i)                           &
[1575]332                      + f1_mg_b(k,l) * p_mg(k,j,i)
333             ENDDO
334          ENDDO
[1898]335       ENDDO
336       !$OMP END PARALLEL
[1575]337!
338!--    Horizontal boundary conditions
339       CALL exchange_horiz( r, 1)
340
341       IF ( .NOT. bc_lr_cyc )  THEN
[3182]342          IF ( bc_dirichlet_l  .OR.  bc_radiation_l )  THEN
[1762]343             r(:,:,nxl_mg(l)-1) = r(:,:,nxl_mg(l))
344          ENDIF
[3182]345          IF ( bc_dirichlet_r  .OR.  bc_radiation_r )  THEN
[1762]346             r(:,:,nxr_mg(l)+1) = r(:,:,nxr_mg(l))
347          ENDIF
[1575]348       ENDIF
349
350       IF ( .NOT. bc_ns_cyc )  THEN
[3182]351          IF ( bc_dirichlet_n  .OR.  bc_radiation_n )  THEN
[1762]352             r(:,nyn_mg(l)+1,:) = r(:,nyn_mg(l),:)
353          ENDIF
[3182]354          IF ( bc_dirichlet_s  .OR.  bc_radiation_s )  THEN
[1762]355             r(:,nys_mg(l)-1,:) = r(:,nys_mg(l),:)
356          ENDIF
[1575]357       ENDIF
358
359!
[1898]360!--    Boundary conditions at bottom and top of the domain. Points may be within
[1575]361!--    buildings, but that doesn't matter.
362       IF ( ibc_p_b == 1 )  THEN
[1898]363!
364!--       equivalent to r(nzb,:,: ) = r(nzb+1,:,:)
365          r(nzb,:,: ) = r(ind_even_odd+1,:,:)
[1575]366       ELSE
367          r(nzb,:,: ) = 0.0_wp
368       ENDIF
369
370       IF ( ibc_p_t == 1 )  THEN
[1898]371!
372!--       equivalent to r(nzt_mg(l)+1,:,: ) = r(nzt_mg(l),:,:)
373          r(nzt_mg(l)+1,:,: ) = r(ind_even_odd,:,:)
[1575]374       ELSE
375          r(nzt_mg(l)+1,:,: ) = 0.0_wp
376       ENDIF
377
378       CALL cpu_log( log_point_s(53), 'resid', 'stop' )
379
[1931]380    END SUBROUTINE resid
[1575]381
382
383!------------------------------------------------------------------------------!
384! Description:
385! ------------
[1682]386!> Interpolates the residual on the next coarser grid with "full weighting"
387!> scheme
[1575]388!------------------------------------------------------------------------------!
[1931]389    SUBROUTINE restrict( f_mg, r )
[1575]390
[1682]391
[1575]392       USE control_parameters,                                                 &
[2696]393           ONLY:  bc_lr_cyc, bc_ns_cyc, ibc_p_b, ibc_p_t
[1575]394
395       USE indices,                                                            &
[1898]396           ONLY:  nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
[1575]397
398       IMPLICIT NONE
399
[1898]400       INTEGER(iwp) ::  i    !< index variable along x on finer grid
401       INTEGER(iwp) ::  ic   !< index variable along x on coarser grid
402       INTEGER(iwp) ::  j    !< index variable along y on finer grid
403       INTEGER(iwp) ::  jc   !< index variable along y on coarser grid
404       INTEGER(iwp) ::  k    !< index variable along z on finer grid
405       INTEGER(iwp) ::  kc   !< index variable along z on coarser grid
406       INTEGER(iwp) ::  l    !< index indicating finer grid level
407       INTEGER(iwp) ::  km1  !< index variable along z dimension (k-1 on finer level)
408       INTEGER(iwp) ::  kp1  !< index variable along z dimension (k+1 on finer level)
[1575]409
410
411       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
412                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
[1898]413                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::       &
414                                         f_mg  !< Residual on coarser grid level
[1575]415
416       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level+1)+1,                         &
417                           nys_mg(grid_level+1)-1:nyn_mg(grid_level+1)+1,      &
[1898]418                           nxl_mg(grid_level+1)-1:nxr_mg(grid_level+1)+1) ::   &
419                                         r !< Residual on finer grid level
[1575]420
421!
422!--    Interpolate the residual
423       l = grid_level
424
425       CALL cpu_log( log_point_s(54), 'restrict', 'start' )
426!
[1898]427!--    No wall treatment
428       !$OMP PARALLEL PRIVATE (i,j,k,ic,jc,kc,km1,kp1)
[2073]429       !$OMP DO SCHEDULE( STATIC )
[1898]430       DO  ic = nxl_mg(l), nxr_mg(l)
431          i = 2*ic
432          DO  jc = nys_mg(l), nyn_mg(l)
[1575]433!
[1898]434!--          Calculation for the first point along k
435             j  = 2*jc
436!
437!--          Calculation for the other points along k
438             !DIR$ IVDEP
439             DO k = ind_even_odd+1, nzt_mg(l+1)    ! Fine grid at this point
440                km1 = k-ind_even_odd-1
[1575]441                kp1 = k-ind_even_odd
[1898]442                kc  = k-ind_even_odd               ! Coarse grid index
[1575]443
[1898]444                f_mg(kc,jc,ic) = 1.0_wp / 64.0_wp * (                      &
445                               8.0_wp * r(k,j,i)                            &
446                             + 4.0_wp * ( r(k,j,i-1)     + r(k,j,i+1)     + &
447                                          r(k,j+1,i)     + r(k,j-1,i)     ) &
448                             + 2.0_wp * ( r(k,j-1,i-1)   + r(k,j+1,i-1)   + &
449                                          r(k,j-1,i+1)   + r(k,j+1,i+1)   ) &
450                             + 4.0_wp * r(km1,j,i)                          &
451                             + 2.0_wp * ( r(km1,j,i-1)   + r(km1,j,i+1)   + &
452                                          r(km1,j+1,i)   + r(km1,j-1,i)   ) &
453                             +          ( r(km1,j-1,i-1) + r(km1,j+1,i-1) + &
454                                          r(km1,j-1,i+1) + r(km1,j+1,i+1) ) &
455                             + 4.0_wp * r(kp1,j,i)                          &
456                             + 2.0_wp * ( r(kp1,j,i-1)   + r(kp1,j,i+1)   + &
457                                          r(kp1,j+1,i)   + r(kp1,j-1,i)   ) &
458                             +          ( r(kp1,j-1,i-1) + r(kp1,j+1,i-1) + &
459                                          r(kp1,j-1,i+1) + r(kp1,j+1,i+1) ) &
460                                        )
[1575]461             ENDDO
462          ENDDO
[1898]463       ENDDO
[2073]464       !$OMP ENDDO
[1898]465       !$OMP END PARALLEL
[1575]466
467!
[1898]468!--    Ghost point exchange
469       CALL exchange_horiz( f_mg, 1)
[1575]470!
471!--    Horizontal boundary conditions
472       IF ( .NOT. bc_lr_cyc )  THEN
[3182]473          IF ( bc_dirichlet_l  .OR.  bc_radiation_l )  THEN
[1762]474             f_mg(:,:,nxl_mg(l)-1) = f_mg(:,:,nxl_mg(l))
475          ENDIF
[3182]476          IF ( bc_dirichlet_r  .OR.  bc_radiation_r )  THEN
[1762]477             f_mg(:,:,nxr_mg(l)+1) = f_mg(:,:,nxr_mg(l))
478          ENDIF
[1575]479       ENDIF
480
481       IF ( .NOT. bc_ns_cyc )  THEN
[3182]482          IF ( bc_dirichlet_n  .OR.  bc_radiation_n )  THEN
[1762]483             f_mg(:,nyn_mg(l)+1,:) = f_mg(:,nyn_mg(l),:)
484          ENDIF
[3182]485          IF ( bc_dirichlet_s  .OR.  bc_radiation_s )  THEN
[1762]486             f_mg(:,nys_mg(l)-1,:) = f_mg(:,nys_mg(l),:)
487          ENDIF
[1575]488       ENDIF
489
490!
491!--    Boundary conditions at bottom and top of the domain.
492!--    These points are not handled by the above loop. Points may be within
[1898]493!--    buildings, but that doesn't matter. Remark: f_mg is ordered sequentielly
494!--    after interpolation on coarse grid (is ordered in odd-even blocks further
495!--    below).
[1575]496       IF ( ibc_p_b == 1 )  THEN
497          f_mg(nzb,:,: ) = f_mg(nzb+1,:,:)
498       ELSE
499          f_mg(nzb,:,: ) = 0.0_wp
500       ENDIF
501
502       IF ( ibc_p_t == 1 )  THEN
503          f_mg(nzt_mg(l)+1,:,: ) = f_mg(nzt_mg(l),:,:)
504       ELSE
505          f_mg(nzt_mg(l)+1,:,: ) = 0.0_wp
506       ENDIF
507
508       CALL cpu_log( log_point_s(54), 'restrict', 'stop' )
509!
[1898]510!--    Since residual is in sequential order after interpolation, an additional
511!--    sorting in odd-even blocks along z dimension is required at this point.
[1575]512       CALL sort_k_to_even_odd_blocks( f_mg , l)
513
[1931]514    END SUBROUTINE restrict
[1575]515
516
517!------------------------------------------------------------------------------!
518! Description:
519! ------------
[1682]520!> Interpolates the correction of the perturbation pressure
521!> to the next finer grid.
[1575]522!------------------------------------------------------------------------------!
[1931]523    SUBROUTINE prolong( p, temp )
[1575]524
[1682]525
[1575]526       USE control_parameters,                                                 &
[2696]527           ONLY:  bc_lr_cyc, bc_ns_cyc, ibc_p_b, ibc_p_t
[1575]528       USE indices,                                                            &
529           ONLY:  nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
530
531       IMPLICIT NONE
532
[1898]533       INTEGER(iwp) ::  i   !< index variable along x on coarser grid level
534       INTEGER(iwp) ::  j   !< index variable along y on coarser grid level
535       INTEGER(iwp) ::  k   !< index variable along z on coarser grid level
536       INTEGER(iwp) ::  l   !< index indicating finer grid level
537       INTEGER(iwp) ::  kp1 !< index variable along z
538       INTEGER(iwp) ::  ke  !< index for prolog even
539       INTEGER(iwp) ::  ko  !< index for prolog odd
[1575]540
541       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1,                         &
542                           nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1,      &
[1898]543                           nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1 ) ::  &
544                               p     !< perturbation pressure on coarser grid level
[1575]545
546       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
547                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
[1898]548                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::       &
549                               temp  !< perturbation pressure on finer grid level
[1575]550
551
552       CALL cpu_log( log_point_s(55), 'prolong', 'start' )
553
554!
555!--    First, store elements of the coarser grid on the next finer grid
556       l = grid_level
557       ind_even_odd = even_odd_level(grid_level-1)
558
559       !$OMP PARALLEL PRIVATE (i,j,k,kp1,ke,ko)
560       !$OMP DO
561       DO  i = nxl_mg(l-1), nxr_mg(l-1)
562          DO  j = nys_mg(l-1), nyn_mg(l-1)
563
564             !DIR$ IVDEP
565             DO k = ind_even_odd+1, nzt_mg(l-1)
566                kp1 = k - ind_even_odd
567                ke  = 2 * ( k-ind_even_odd - 1 ) + 1
568                ko  = 2 * k - 1
569!
570!--             Points of the coarse grid are directly stored on the next finer
571!--             grid
572                temp(ko,2*j,2*i)   = p(k,j,i)
573!
574!--             Points between two coarse-grid points
575                temp(ko,2*j,2*i+1) = 0.5_wp * ( p(k,j,i) + p(k,j,i+1) )
576                temp(ko,2*j+1,2*i) = 0.5_wp * ( p(k,j,i) + p(k,j+1,i) )
577                temp(ke,2*j,2*i)   = 0.5_wp * ( p(k,j,i) + p(kp1,j,i) )
578!
579!--             Points in the center of the planes stretched by four points
580!--             of the coarse grid cube
581                temp(ko,2*j+1,2*i+1) = 0.25_wp * ( p(k,j,i)   + p(k,j,i+1) +   &
582                                                   p(k,j+1,i) + p(k,j+1,i+1) )
583                temp(ke,2*j,2*i+1)   = 0.25_wp * ( p(k,j,i)   + p(k,j,i+1) +   &
584                                                   p(kp1,j,i) + p(kp1,j,i+1) )
585                temp(ke,2*j+1,2*i)   = 0.25_wp * ( p(k,j,i)   + p(k,j+1,i) +   &
586                                                   p(kp1,j,i) + p(kp1,j+1,i) )
587!
588!--             Points in the middle of coarse grid cube
589                temp(ke,2*j+1,2*i+1) = 0.125_wp *                              &
590                                               ( p(k,j,i)     + p(k,j,i+1)   + &
591                                                 p(k,j+1,i)   + p(k,j+1,i+1) + &
592                                                 p(kp1,j,i)   + p(kp1,j,i+1) + &
593                                                 p(kp1,j+1,i) + p(kp1,j+1,i+1) )
[1898]594
[1575]595             ENDDO
596
597             !DIR$ IVDEP
598             DO k = nzb+1, ind_even_odd
599                kp1 = k + ind_even_odd + 1
600                ke  = 2 * k
601                ko  = 2 * ( k + ind_even_odd )
602!
603!--             Points of the coarse grid are directly stored on the next finer
604!--             grid
605                temp(ko,2*j,2*i)   = p(k,j,i)
606!
607!--             Points between two coarse-grid points
608                temp(ko,2*j,2*i+1) = 0.5_wp * ( p(k,j,i) + p(k,j,i+1) )
609                temp(ko,2*j+1,2*i) = 0.5_wp * ( p(k,j,i) + p(k,j+1,i) )
610                temp(ke,2*j,2*i)   = 0.5_wp * ( p(k,j,i) + p(kp1,j,i) )
611!
612!--             Points in the center of the planes stretched by four points
613!--             of the coarse grid cube
614                temp(ko,2*j+1,2*i+1) = 0.25_wp * ( p(k,j,i)   + p(k,j,i+1) +   &
615                                                   p(k,j+1,i) + p(k,j+1,i+1) )
616                temp(ke,2*j,2*i+1)   = 0.25_wp * ( p(k,j,i)   + p(k,j,i+1) +   &
617                                                   p(kp1,j,i) + p(kp1,j,i+1) )
618                temp(ke,2*j+1,2*i)   = 0.25_wp * ( p(k,j,i)   + p(k,j+1,i) +   &
619                                                   p(kp1,j,i) + p(kp1,j+1,i) )
620!
621!--             Points in the middle of coarse grid cube
622                temp(ke,2*j+1,2*i+1) = 0.125_wp *                              &
623                                               ( p(k,j,i)     + p(k,j,i+1)   + &
624                                                 p(k,j+1,i)   + p(k,j+1,i+1) + &
625                                                 p(kp1,j,i)   + p(kp1,j,i+1) + &
626                                                 p(kp1,j+1,i) + p(kp1,j+1,i+1) )
[1898]627
[1575]628             ENDDO
629
630          ENDDO
631       ENDDO
632       !$OMP END PARALLEL
633
634       ind_even_odd = even_odd_level(grid_level)
635!
636!--    Horizontal boundary conditions
637       CALL exchange_horiz( temp, 1)
638
639       IF ( .NOT. bc_lr_cyc )  THEN
[3182]640          IF ( bc_dirichlet_l  .OR.  bc_radiation_l )  THEN
[1762]641             temp(:,:,nxl_mg(l)-1) = temp(:,:,nxl_mg(l))
642          ENDIF
[3182]643          IF ( bc_dirichlet_r  .OR.  bc_radiation_r )  THEN
[1762]644             temp(:,:,nxr_mg(l)+1) = temp(:,:,nxr_mg(l))
645          ENDIF
[1575]646       ENDIF
647
648       IF ( .NOT. bc_ns_cyc )  THEN
[3182]649          IF ( bc_dirichlet_n  .OR.  bc_radiation_n )  THEN
[1762]650             temp(:,nyn_mg(l)+1,:) = temp(:,nyn_mg(l),:)
651          ENDIF
[3182]652          IF ( bc_dirichlet_s  .OR.  bc_radiation_s )  THEN
[1762]653             temp(:,nys_mg(l)-1,:) = temp(:,nys_mg(l),:)
654          ENDIF
[1575]655       ENDIF
656
657!
658!--    Bottom and top boundary conditions
659       IF ( ibc_p_b == 1 )  THEN
[1898]660!
661!--       equivalent to temp(nzb,:,: ) = temp(nzb+1,:,:)
[1575]662          temp(nzb,:,: ) = temp(ind_even_odd+1,:,:)
663       ELSE
664          temp(nzb,:,: ) = 0.0_wp
665       ENDIF
666
667       IF ( ibc_p_t == 1 )  THEN
[1898]668!
669!--       equivalent to temp(nzt_mg(l)+1,:,: ) = temp(nzt_mg(l),:,:)
[1575]670          temp(nzt_mg(l)+1,:,: ) = temp(ind_even_odd,:,:)
671       ELSE
672          temp(nzt_mg(l)+1,:,: ) = 0.0_wp
673       ENDIF
674
675       CALL cpu_log( log_point_s(55), 'prolong', 'stop' )
676
[1931]677    END SUBROUTINE prolong
[1575]678
679
680!------------------------------------------------------------------------------!
681! Description:
682! ------------
[1682]683!> Relaxation method for the multigrid scheme. A Gauss-Seidel iteration with
684!> 3D-Red-Black decomposition (GS-RB) is used.
[1575]685!------------------------------------------------------------------------------!
[1931]686    SUBROUTINE redblack( f_mg, p_mg )
[1575]687
[1682]688
[1575]689       USE arrays_3d,                                                          &
[3241]690           ONLY:  rho_air_mg
[1575]691
692       USE control_parameters,                                                 &
[2696]693           ONLY:  bc_lr_cyc, bc_ns_cyc, ibc_p_b, ibc_p_t, ngsrb 
[1575]694
695       USE grid_variables,                                                     &
696           ONLY:  ddx2_mg, ddy2_mg
697
698       USE indices,                                                            &
[1898]699           ONLY:  nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
[1575]700
701       IMPLICIT NONE
702
[1898]703       INTEGER(iwp) :: color    !< grid point color, either red or black
704       INTEGER(iwp) :: i        !< index variable along x
705       INTEGER(iwp) :: ic       !< index variable along x
706       INTEGER(iwp) :: j        !< index variable along y
707       INTEGER(iwp) :: jc       !< index variable along y
708       INTEGER(iwp) :: jj       !< index variable along y
709       INTEGER(iwp) :: k        !< index variable along z
710       INTEGER(iwp) :: l        !< grid level
711       INTEGER(iwp) :: n        !< loop variable Gauß-Seidel iterations
712       INTEGER(iwp) :: km1      !< index variable (k-1)
713       INTEGER(iwp) :: kp1      !< index variable (k+1)
[1575]714
[1898]715       LOGICAL      :: unroll   !< flag indicating whether loop unrolling is possible
[1575]716
717       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
718                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
[1898]719                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::       &
720                                      f_mg  !< residual of perturbation pressure
[1575]721       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
722                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
[1898]723                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::       &
724                                      p_mg  !< perturbation pressure
[1575]725
726       l = grid_level
727
[1898]728       unroll = ( MOD( nyn_mg(l)-nys_mg(l)+1, 4 ) == 0  .AND.                  &
[1575]729                  MOD( nxr_mg(l)-nxl_mg(l)+1, 2 ) == 0 )
730
731       DO  n = 1, ngsrb
732       
733          DO  color = 1, 2
734
[1898]735             IF ( .NOT. unroll )  THEN
[1575]736
[1898]737                CALL cpu_log( log_point_s(36), 'redblack_no_unroll_f', 'start' )
[1575]738!
[1898]739!--             Without unrolling of loops, no cache optimization
740                !$OMP PARALLEL PRIVATE (i,j,k,km1,kp1)
741                !$OMP DO
742                DO  i = nxl_mg(l), nxr_mg(l), 2
743                   DO  j = nys_mg(l) + 2 - color, nyn_mg(l), 2
744                      !DIR$ IVDEP
745                      DO  k = ind_even_odd+1, nzt_mg(l)
746                         km1 = k-ind_even_odd-1
747                         kp1 = k-ind_even_odd
748                         p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
[2232]749                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
[1898]750                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
[2232]751                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
[1898]752                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
[1575]753                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
754                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
[1898]755                               - f_mg(k,j,i)                   )
[1575]756                      ENDDO
757                   ENDDO
[1898]758                ENDDO
[1575]759   
[1898]760                !$OMP DO
761                DO  i = nxl_mg(l)+1, nxr_mg(l), 2
762                   DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
763                       !DIR$ IVDEP
764                       DO  k = ind_even_odd+1, nzt_mg(l)
765                         km1 = k-ind_even_odd-1
766                         kp1 = k-ind_even_odd
767                         p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
[2232]768                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
[1898]769                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
[2232]770                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
[1575]771                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
772                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
773                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
[1898]774                               - f_mg(k,j,i)                   )
[1575]775                      ENDDO
776                   ENDDO
[1898]777                ENDDO
[1575]778 
[1898]779                !$OMP DO
780                DO  i = nxl_mg(l), nxr_mg(l), 2
781                   DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
782                      !DIR$ IVDEP
783                      DO  k = nzb+1, ind_even_odd
784                         km1 = k+ind_even_odd
785                         kp1 = k+ind_even_odd+1
786                         p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
[2232]787                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
[1575]788                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
[2232]789                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
[1575]790                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
791                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
792                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
[1898]793                               - f_mg(k,j,i)                   )
[1575]794                      ENDDO
795                   ENDDO
[1898]796                ENDDO
[1575]797
[1898]798                !$OMP DO
799                DO  i = nxl_mg(l)+1, nxr_mg(l), 2
800                   DO  j = nys_mg(l) + 2 - color, nyn_mg(l), 2
801                      !DIR$ IVDEP
802                      DO  k = nzb+1, ind_even_odd
803                         km1 = k+ind_even_odd
804                         kp1 = k+ind_even_odd+1
805                         p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
[2232]806                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
[1575]807                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
[2232]808                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
[1575]809                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
810                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
811                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
[1898]812                               - f_mg(k,j,i)                   )
[1575]813                      ENDDO
814                   ENDDO
[1898]815                ENDDO
816                !$OMP END PARALLEL
[1575]817
[1898]818                CALL cpu_log( log_point_s(36), 'redblack_no_unroll_f', 'stop' )
[1575]819
[1898]820             ELSE
[1575]821!
[1898]822!--              Loop unrolling along y, only one i loop for better cache use
823                CALL cpu_log( log_point_s(38), 'redblack_unroll_f', 'start' )
[1575]824
[1898]825                !$OMP PARALLEL PRIVATE (i,j,k,ic,jc,km1,kp1,jj)
826                !$OMP DO
827                DO  ic = nxl_mg(l), nxr_mg(l), 2
828                   DO  jc = nys_mg(l), nyn_mg(l), 4
829                      i  = ic
830                      jj = jc+2-color
831                      !DIR$ IVDEP
832                      DO  k = ind_even_odd+1, nzt_mg(l)
833                         km1 = k-ind_even_odd-1
834                         kp1 = k-ind_even_odd
835                         j   = jj
836                         p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
[2232]837                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
[1575]838                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
[2232]839                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
[1575]840                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
841                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
842                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
[1898]843                               - f_mg(k,j,i)                   )
844                         j = jj+2
845                         p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
[2232]846                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
[1575]847                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
[2232]848                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
[1575]849                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
850                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
851                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
[1898]852                               - f_mg(k,j,i)                   )
853                      ENDDO
[1575]854
[1898]855                      i  = ic+1
856                      jj = jc+color-1
857                      !DIR$ IVDEP
858                      DO  k = ind_even_odd+1, nzt_mg(l)
859                         km1 = k-ind_even_odd-1
860                         kp1 = k-ind_even_odd
861                         j   = jj
862                         p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
[2232]863                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
[1898]864                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
[2232]865                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
[1575]866                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
867                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
868                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
[1898]869                               - f_mg(k,j,i)                   )
870                         j = jj+2
871                         p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
[2232]872                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
[1898]873                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
[2232]874                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
[1898]875                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
[1575]876                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
877                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
[1898]878                               - f_mg(k,j,i)                   )
879                      ENDDO
[1575]880
[1898]881                      i  = ic
882                      jj = jc+color-1
883                      !DIR$ IVDEP
884                      DO  k = nzb+1, ind_even_odd
885                         km1 = k+ind_even_odd
886                         kp1 = k+ind_even_odd+1
887                         j   = jj
888                         p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
[2232]889                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
[1898]890                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
[2232]891                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
[1898]892                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
[1575]893                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
894                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
[1898]895                               - f_mg(k,j,i)                   )
896                         j = jj+2
897                         p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
[2232]898                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
[1575]899                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
[2232]900                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
[1575]901                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
902                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
903                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
[1898]904                               - f_mg(k,j,i)                   )
905                      ENDDO
[1575]906
[1898]907                      i  = ic+1
908                      jj = jc+2-color
909                      !DIR$ IVDEP
910                      DO  k = nzb+1, ind_even_odd
911                         km1 = k+ind_even_odd
912                         kp1 = k+ind_even_odd+1
913                         j   = jj
914                         p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
[2232]915                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
[1898]916                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
[2232]917                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
[1575]918                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
919                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
920                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
[1898]921                               - f_mg(k,j,i)                   )
922                         j = jj+2
923                         p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
[2232]924                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
[1575]925                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
[2232]926                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
[1898]927                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
[1575]928                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
929                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
[1898]930                               - f_mg(k,j,i)                   )
[1575]931                      ENDDO
932
933                   ENDDO
[1898]934                ENDDO
935                !$OMP END PARALLEL
[1575]936
[1898]937                CALL cpu_log( log_point_s(38), 'redblack_unroll_f', 'stop' )
[1575]938
939             ENDIF
940
941!
942!--          Horizontal boundary conditions
943             CALL special_exchange_horiz( p_mg, color )
944
[1762]945             IF ( .NOT. bc_lr_cyc )  THEN
[3182]946                IF ( bc_dirichlet_l  .OR.  bc_radiation_l )  THEN
[1575]947                   p_mg(:,:,nxl_mg(l)-1) = p_mg(:,:,nxl_mg(l))
948                ENDIF
[3182]949                IF ( bc_dirichlet_r  .OR.  bc_radiation_r )  THEN
[1575]950                   p_mg(:,:,nxr_mg(l)+1) = p_mg(:,:,nxr_mg(l))
951                ENDIF
952             ENDIF
953
954             IF ( .NOT. bc_ns_cyc )  THEN
[3182]955                IF ( bc_dirichlet_n  .OR.  bc_radiation_n )  THEN
[1575]956                   p_mg(:,nyn_mg(l)+1,:) = p_mg(:,nyn_mg(l),:)
957                ENDIF
[3182]958                IF ( bc_dirichlet_s  .OR.  bc_radiation_s )  THEN
[1575]959                   p_mg(:,nys_mg(l)-1,:) = p_mg(:,nys_mg(l),:)
960                ENDIF
961             ENDIF
962
963!
964!--          Bottom and top boundary conditions
965             IF ( ibc_p_b == 1 )  THEN
[1898]966!
967!--             equivalent to p_mg(nzb,:,: ) = p_mg(nzb+1,:,:)
[1575]968                p_mg(nzb,:,: ) = p_mg(ind_even_odd+1,:,:)
969             ELSE
970                p_mg(nzb,:,: ) = 0.0_wp
971             ENDIF
972
973             IF ( ibc_p_t == 1 )  THEN
[1898]974!
975!--             equivalent to p_mg(nzt_mg(l)+1,:,: ) = p_mg(nzt_mg(l),:,:)
[1575]976                p_mg(nzt_mg(l)+1,:,: ) = p_mg(ind_even_odd,:,:)
977             ELSE
978                p_mg(nzt_mg(l)+1,:,: ) = 0.0_wp
979             ENDIF
980
981          ENDDO
[1898]982
[1575]983       ENDDO
984
[1931]985    END SUBROUTINE redblack
[1575]986
987
988!------------------------------------------------------------------------------!
989! Description:
990! ------------
[1682]991!> Sort k-Dimension from sequential into blocks of even and odd.
992!> This is required to vectorize the red-black subroutine.
993!> Version for 3D-REAL arrays
[1575]994!------------------------------------------------------------------------------!
[1682]995    SUBROUTINE sort_k_to_even_odd_blocks( p_mg , glevel )
[1575]996
[1682]997
[1575]998       USE indices,                                                            &
999           ONLY:  nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
1000
1001       IMPLICIT NONE
1002
[1898]1003       INTEGER(iwp), INTENT(IN) ::  glevel  !< grid level
[1575]1004
1005       REAL(wp), DIMENSION(nzb:nzt_mg(glevel)+1,                               &
1006                           nys_mg(glevel)-1:nyn_mg(glevel)+1,                  &
[1898]1007                           nxl_mg(glevel)-1:nxr_mg(glevel)+1) ::               &
1008                                      p_mg  !< array to be sorted
[1575]1009!
1010!--    Local variables
[1898]1011       INTEGER(iwp) :: i        !< index variable along x
1012       INTEGER(iwp) :: j        !< index variable along y
1013       INTEGER(iwp) :: k        !< index variable along z
1014       INTEGER(iwp) :: l        !< grid level
1015       INTEGER(iwp) :: ind      !< index variable along z
1016       REAL(wp), DIMENSION(nzb:nzt_mg(glevel)+1) ::  tmp  !< odd-even sorted temporary array
[1575]1017
1018
1019       CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'start' )
1020
1021       l = glevel
1022       ind_even_odd = even_odd_level(l)
1023
1024       !$OMP PARALLEL PRIVATE (i,j,k,ind,tmp)
1025       !$OMP DO
1026       DO  i = nxl_mg(l)-1, nxr_mg(l)+1
1027          DO  j = nys_mg(l)-1, nyn_mg(l)+1
1028
1029!
1030!--          Sort the data with even k index
1031             ind = nzb-1
1032             DO  k = nzb, nzt_mg(l), 2
1033                ind = ind + 1
1034                tmp(ind) = p_mg(k,j,i)
1035             ENDDO
1036!
1037!--          Sort the data with odd k index
1038             DO  k = nzb+1, nzt_mg(l)+1, 2
1039                ind = ind + 1
1040                tmp(ind) = p_mg(k,j,i)
1041             ENDDO
1042
1043             p_mg(:,j,i) = tmp
1044
1045          ENDDO
1046       ENDDO
1047       !$OMP END PARALLEL
1048
1049       CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'stop' )
1050
1051    END SUBROUTINE sort_k_to_even_odd_blocks
1052
1053
1054!------------------------------------------------------------------------------!
1055! Description:
1056! ------------
[1682]1057!> Sort k-Dimension from sequential into blocks of even and odd.
1058!> This is required to vectorize the red-black subroutine.
1059!> Version for 1D-REAL arrays
[1575]1060!------------------------------------------------------------------------------!
[1682]1061    SUBROUTINE sort_k_to_even_odd_blocks_1d( f_mg, f_mg_b, glevel )
[1575]1062
[1682]1063
[1575]1064       USE indices,                                                            &
1065           ONLY:  nzb, nzt_mg
1066
1067       IMPLICIT NONE
1068
[1898]1069       INTEGER(iwp), INTENT(IN) ::  glevel  !< grid level
[1575]1070
[1898]1071       REAL(wp), DIMENSION(nzb+1:nzt_mg(glevel)) ::  f_mg    !< 1D input array
1072       REAL(wp), DIMENSION(nzb:nzt_mg(glevel)+1) ::  f_mg_b  !< 1D output array
[1575]1073
1074!
1075!--    Local variables
[1898]1076       INTEGER(iwp) :: ind   !< index variable along z
1077       INTEGER(iwp) :: k     !< index variable along z
[1575]1078
1079
1080       ind = nzb - 1
1081!
1082!--    Sort the data with even k index
1083       DO  k = nzb, nzt_mg(glevel), 2
1084          ind = ind + 1
1085          IF ( k >= nzb+1  .AND.  k <= nzt_mg(glevel) )  THEN
1086             f_mg_b(ind) = f_mg(k)
1087          ENDIF
1088       ENDDO
1089!
1090!--    Sort the data with odd k index
1091       DO  k = nzb+1, nzt_mg(glevel)+1, 2
1092          ind = ind + 1
1093          IF( k >= nzb+1  .AND.  k <= nzt_mg(glevel) )  THEN
1094             f_mg_b(ind) = f_mg(k)
1095          ENDIF
1096       ENDDO
1097
1098    END SUBROUTINE sort_k_to_even_odd_blocks_1d
1099
1100
1101!------------------------------------------------------------------------------!
1102! Description:
1103! ------------
[1682]1104!> Sort k-Dimension from sequential into blocks of even and odd.
1105!> This is required to vectorize the red-black subroutine.
1106!> Version for 2D-INTEGER arrays
[1575]1107!------------------------------------------------------------------------------!
[3725]1108!    SUBROUTINE sort_k_to_even_odd_blocks_int( i_mg , glevel )
[1575]1109!
[3725]1110!
1111!       USE indices,                                                            &
1112!           ONLY:  nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
1113!
1114!       IMPLICIT NONE
1115!
1116!       INTEGER(iwp), INTENT(IN) ::  glevel  !< grid level
1117!
1118!       INTEGER(iwp), DIMENSION(nzb:nzt_mg(glevel)+1,                           &
1119!                               nys_mg(glevel)-1:nyn_mg(glevel)+1,              &
1120!                               nxl_mg(glevel)-1:nxr_mg(glevel)+1) ::           &
1121!                                    i_mg    !< array to be sorted
1122!
[1575]1123!--    Local variables
[3725]1124!       INTEGER(iwp) :: i        !< index variabel along x
1125!       INTEGER(iwp) :: j        !< index variable along y
1126!       INTEGER(iwp) :: k        !< index variable along z
1127!       INTEGER(iwp) :: l        !< grid level
1128!       INTEGER(iwp) :: ind      !< index variable along z
1129!       INTEGER(iwp),DIMENSION(nzb:nzt_mg(glevel)+1) ::  tmp  !< temporary odd-even sorted array
[1575]1130!
[3725]1131!
1132!       CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'start' )
1133!
1134!       l = glevel
1135!       ind_even_odd = even_odd_level(l)
1136!
1137!       DO  i = nxl_mg(l)-1, nxr_mg(l)+1
1138!          DO  j = nys_mg(l)-1, nyn_mg(l)+1
1139!
1140!
[1575]1141!--          Sort the data with even k index
[3725]1142!             ind = nzb-1
1143!             DO  k = nzb, nzt_mg(l), 2
1144!                ind = ind + 1
1145!                tmp(ind) = i_mg(k,j,i)
1146!             ENDDO
[1575]1147!
1148!--          Sort the data with odd k index
[3725]1149!             DO  k = nzb+1, nzt_mg(l)+1, 2
1150!                ind = ind + 1
1151!                tmp(ind) = i_mg(k,j,i)
1152!             ENDDO
1153!
1154!             i_mg(:,j,i) = tmp
1155!
1156!          ENDDO
1157!       ENDDO
1158!
1159!       CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'stop' )
1160!
1161!    END SUBROUTINE sort_k_to_even_odd_blocks_int
[1575]1162
1163
1164!------------------------------------------------------------------------------!
1165! Description:
1166! ------------
[1682]1167!> Sort k-dimension from blocks of even and odd into sequential
[1575]1168!------------------------------------------------------------------------------!
[1682]1169    SUBROUTINE sort_k_to_sequential( p_mg )
[1575]1170
[1682]1171
[1575]1172       USE control_parameters,                                                 &
1173           ONLY:  grid_level
1174
1175       USE indices,                                                            &
1176           ONLY:  nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
1177
1178       IMPLICIT NONE
1179
1180       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
1181                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
[1898]1182                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::       &
1183                                     p_mg  !< array to be sorted
[1575]1184!
1185!--    Local variables
[1898]1186       INTEGER(iwp) :: i        !< index variable along x
1187       INTEGER(iwp) :: j        !< index variable along y
1188       INTEGER(iwp) :: k        !< index variable along z
1189       INTEGER(iwp) :: l        !< grid level
1190       INTEGER(iwp) :: ind      !< index variable along z
[1575]1191
1192       REAL(wp),DIMENSION(nzb:nzt_mg(grid_level)+1) ::  tmp
1193
1194
1195       l = grid_level
1196
1197       !$OMP PARALLEL PRIVATE (i,j,k,ind,tmp)
1198       !$OMP DO
1199       DO  i = nxl_mg(l)-1, nxr_mg(l)+1
1200          DO  j = nys_mg(l)-1, nyn_mg(l)+1
1201
1202             ind = nzb - 1
1203             tmp = p_mg(:,j,i)
1204             DO  k = nzb, nzt_mg(l), 2
1205                ind = ind + 1
1206                p_mg(k,j,i) = tmp(ind)
1207             ENDDO
1208
1209             DO  k = nzb+1, nzt_mg(l)+1, 2
1210                ind = ind + 1
1211                p_mg(k,j,i) = tmp(ind)
1212             ENDDO
1213          ENDDO
1214       ENDDO
1215       !$OMP END PARALLEL
1216
1217    END SUBROUTINE sort_k_to_sequential
1218
1219
[1682]1220!------------------------------------------------------------------------------!
1221! Description:
1222! ------------
1223!> Gather subdomain data from all PEs.
1224!------------------------------------------------------------------------------!
[4429]1225#if defined( __parallel )
[1931]1226    SUBROUTINE mg_gather( f2, f2_sub )
[1575]1227
1228       USE control_parameters,                                                 &
1229           ONLY:  grid_level
1230
1231       USE cpulog,                                                             &
1232           ONLY:  cpu_log, log_point_s
1233
1234       USE indices,                                                            &
1235           ONLY:  mg_loc_ind, nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
1236
1237       IMPLICIT NONE
1238
[1682]1239       INTEGER(iwp) ::  i       !<
1240       INTEGER(iwp) ::  il      !<
1241       INTEGER(iwp) ::  ir      !<
1242       INTEGER(iwp) ::  j       !<
1243       INTEGER(iwp) ::  jn      !<
1244       INTEGER(iwp) ::  js      !<
1245       INTEGER(iwp) ::  k       !<
1246       INTEGER(iwp) ::  nwords  !<
[1575]1247
1248       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
1249                       nys_mg(grid_level)-1:nyn_mg(grid_level)+1,              &
[1682]1250                       nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f2    !<
[1575]1251       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
1252                       nys_mg(grid_level)-1:nyn_mg(grid_level)+1,              &
[1682]1253                       nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f2_l  !<
[1575]1254
1255       REAL(wp), DIMENSION(nzb:mg_loc_ind(5,myid)+1,                           &
1256                           mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1,          &
[1682]1257                           mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) ::  f2_sub  !<
[1575]1258
1259
1260       CALL cpu_log( log_point_s(34), 'mg_gather', 'start' )
1261
1262       f2_l = 0.0_wp
1263
1264!
1265!--    Store the local subdomain array on the total array
1266       js = mg_loc_ind(3,myid)
1267       IF ( south_border_pe )  js = js - 1
1268       jn = mg_loc_ind(4,myid)
1269       IF ( north_border_pe )  jn = jn + 1
1270       il = mg_loc_ind(1,myid)
1271       IF ( left_border_pe )   il = il - 1
1272       ir = mg_loc_ind(2,myid)
1273       IF ( right_border_pe )  ir = ir + 1
1274       DO  i = il, ir
1275          DO  j = js, jn
1276             DO  k = nzb, nzt_mg(grid_level)+1
1277                f2_l(k,j,i) = f2_sub(k,j,i)
1278             ENDDO
1279          ENDDO
1280       ENDDO
1281
1282!
1283!--    Find out the number of array elements of the total array
1284       nwords = SIZE( f2 )
1285
1286!
1287!--    Gather subdomain data from all PEs
1288       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1289       CALL MPI_ALLREDUCE( f2_l(nzb,nys_mg(grid_level)-1,nxl_mg(grid_level)-1), &
1290                           f2(nzb,nys_mg(grid_level)-1,nxl_mg(grid_level)-1),   &
1291                           nwords, MPI_REAL, MPI_SUM, comm2d, ierr )
1292
1293       CALL cpu_log( log_point_s(34), 'mg_gather', 'stop' )
1294   
[1931]1295    END SUBROUTINE mg_gather
[4429]1296#endif
[1575]1297
1298
[1682]1299!------------------------------------------------------------------------------!
1300! Description:
1301! ------------
1302!> @todo It might be possible to improve the speed of this routine by using
1303!>       non-blocking communication
1304!------------------------------------------------------------------------------!
[4429]1305#if defined( __parallel )
[1931]1306    SUBROUTINE mg_scatter( p2, p2_sub )
[1575]1307
1308       USE control_parameters,                                                 &
1309           ONLY:  grid_level
1310
1311       USE cpulog,                                                             &
1312           ONLY:  cpu_log, log_point_s
1313
1314       USE indices,                                                            &
1315           ONLY:  mg_loc_ind, nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
1316
1317       IMPLICIT NONE
1318
1319       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1,                         &
1320                           nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1,      &
[1682]1321                           nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) ::  p2  !<
[1575]1322
1323       REAL(wp), DIMENSION(nzb:mg_loc_ind(5,myid)+1,                           &
1324                           mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1,          &
[1682]1325                           mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) ::  p2_sub  !<
[1575]1326
1327
1328       CALL cpu_log( log_point_s(35), 'mg_scatter', 'start' )
1329
1330       p2_sub = p2(:,mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, &
1331                     mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1)
1332
1333       CALL cpu_log( log_point_s(35), 'mg_scatter', 'stop' )
1334   
[1931]1335    END SUBROUTINE mg_scatter
[4429]1336#endif
[1575]1337
1338!------------------------------------------------------------------------------!
1339! Description:
1340! ------------
[1682]1341!> This is where the multigrid technique takes place. V- and W- Cycle are
1342!> implemented and steered by the parameter "gamma". Parameter "nue" determines
1343!> the convergence of the multigrid iterative solution. There are nue times
1344!> RB-GS iterations. It should be set to "1" or "2", considering the time effort
1345!> one would like to invest. Last choice shows a very good converging factor,
1346!> but leads to an increase in computing time.
[1575]1347!------------------------------------------------------------------------------!
[1931]1348    RECURSIVE SUBROUTINE next_mg_level( f_mg, p_mg, p3, r )
[1575]1349
1350       USE control_parameters,                                                 &
1351           ONLY:  bc_lr_dirrad, bc_lr_raddir, bc_ns_dirrad, bc_ns_raddir,      &
[3241]1352                  child_domain, gamma_mg, grid_level_count, maximum_grid_level,&
1353                  mg_switch_to_pe0_level, mg_switch_to_pe0, ngsrb
[1575]1354
1355       USE indices,                                                            &
1356           ONLY:  mg_loc_ind, nxl, nxl_mg, nxr, nxr_mg, nys, nys_mg, nyn,      &
1357                  nyn_mg, nzb, nzt, nzt_mg
1358
1359       IMPLICIT NONE
1360
[1898]1361       INTEGER(iwp) ::  i            !< index variable along x
1362       INTEGER(iwp) ::  j            !< index variable along y
1363       INTEGER(iwp) ::  k            !< index variable along z
[1682]1364       INTEGER(iwp) ::  nxl_mg_save  !<
1365       INTEGER(iwp) ::  nxr_mg_save  !<
1366       INTEGER(iwp) ::  nyn_mg_save  !<
1367       INTEGER(iwp) ::  nys_mg_save  !<
1368       INTEGER(iwp) ::  nzt_mg_save  !<
[1575]1369
1370       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
1371                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
[1682]1372                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg  !<
[1575]1373       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
1374                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
[1682]1375                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p_mg  !<
[1575]1376       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
1377                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
[1682]1378                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p3    !<
[1575]1379       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
1380                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
[1682]1381                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: r     !<
[1575]1382
1383       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1,                         &
1384                           nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1,      &
[1682]1385                           nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) ::  f2  !<
[1575]1386       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1,                         &
1387                           nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1,      &
[1682]1388                           nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) ::  p2  !<
[1575]1389
[1682]1390       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  f2_sub  !<
[4429]1391
1392#if defined( __parallel )
[1682]1393       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  p2_sub  !<
[4429]1394#endif
[1575]1395
1396!
1397!--    Restriction to the coarsest grid
1398    10 IF ( grid_level == 1 )  THEN
1399
1400!
1401!--       Solution on the coarsest grid. Double the number of Gauss-Seidel
1402!--       iterations in order to get a more accurate solution.
1403          ngsrb = 2 * ngsrb
1404
1405          ind_even_odd = even_odd_level(grid_level)
1406
[1931]1407          CALL redblack( f_mg, p_mg )
[1575]1408
1409          ngsrb = ngsrb / 2
1410
1411
1412       ELSEIF ( grid_level /= 1 )  THEN
1413
1414          grid_level_count(grid_level) = grid_level_count(grid_level) + 1
1415
1416!
1417!--       Solution on the actual grid level
1418          ind_even_odd = even_odd_level(grid_level)
1419
[1931]1420          CALL redblack( f_mg, p_mg )
[1575]1421
1422!
1423!--       Determination of the actual residual
[1931]1424          CALL resid( f_mg, p_mg, r )
[1575]1425
1426!--       Restriction of the residual (finer grid values!) to the next coarser
1427!--       grid. Therefore, the grid level has to be decremented now. nxl..nzt have
1428!--       to be set to the coarse grid values, because these variables are needed
1429!--       for the exchange of ghost points in routine exchange_horiz
1430          grid_level = grid_level - 1
1431
1432          nxl = nxl_mg(grid_level)
1433          nys = nys_mg(grid_level)
1434          nxr = nxr_mg(grid_level)
1435          nyn = nyn_mg(grid_level)
1436          nzt = nzt_mg(grid_level)
1437
1438          IF ( grid_level == mg_switch_to_pe0_level )  THEN
1439
1440!
1441!--          From this level on, calculations are done on PE0 only.
1442!--          First, carry out restriction on the subdomain.
1443!--          Therefore, indices of the level have to be changed to subdomain
1444!--          values in between (otherwise, the restrict routine would expect
1445!--          the gathered array)
1446
1447             nxl_mg_save = nxl_mg(grid_level)
1448             nxr_mg_save = nxr_mg(grid_level)
1449             nys_mg_save = nys_mg(grid_level)
1450             nyn_mg_save = nyn_mg(grid_level)
1451             nzt_mg_save = nzt_mg(grid_level)
1452             nxl_mg(grid_level) = mg_loc_ind(1,myid)
1453             nxr_mg(grid_level) = mg_loc_ind(2,myid)
1454             nys_mg(grid_level) = mg_loc_ind(3,myid)
1455             nyn_mg(grid_level) = mg_loc_ind(4,myid)
1456             nzt_mg(grid_level) = mg_loc_ind(5,myid)
1457             nxl = mg_loc_ind(1,myid)
1458             nxr = mg_loc_ind(2,myid)
1459             nys = mg_loc_ind(3,myid)
1460             nyn = mg_loc_ind(4,myid)
1461             nzt = mg_loc_ind(5,myid)
1462
1463             ALLOCATE( f2_sub(nzb:nzt_mg(grid_level)+1,                    &
1464                              nys_mg(grid_level)-1:nyn_mg(grid_level)+1,   &
1465                              nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) )
1466
[1931]1467             CALL restrict( f2_sub, r )
[1575]1468
1469!
1470!--          Restore the correct indices of this level
1471             nxl_mg(grid_level) = nxl_mg_save
1472             nxr_mg(grid_level) = nxr_mg_save
1473             nys_mg(grid_level) = nys_mg_save
1474             nyn_mg(grid_level) = nyn_mg_save
1475             nzt_mg(grid_level) = nzt_mg_save
1476             nxl = nxl_mg(grid_level)
1477             nxr = nxr_mg(grid_level)
1478             nys = nys_mg(grid_level)
1479             nyn = nyn_mg(grid_level)
1480             nzt = nzt_mg(grid_level)
1481!
1482!--          Gather all arrays from the subdomains on PE0
[4429]1483#if defined( __parallel )
[1931]1484             CALL mg_gather( f2, f2_sub )
[4429]1485#endif
[1575]1486
1487!
1488!--          Set switch for routine exchange_horiz, that no ghostpoint exchange
1489!--          has to be carried out from now on
1490             mg_switch_to_pe0 = .TRUE.
1491
1492!
1493!--          In case of non-cyclic lateral boundary conditions, both in- and
1494!--          outflow conditions have to be used on all PEs after the switch,
1495!--          because then they have the total domain.
1496             IF ( bc_lr_dirrad )  THEN
[3182]1497                bc_dirichlet_l  = .TRUE.
1498                bc_dirichlet_r  = .FALSE.
1499                bc_radiation_l = .FALSE.
1500                bc_radiation_r = .TRUE.
[1575]1501             ELSEIF ( bc_lr_raddir )  THEN
[3182]1502                bc_dirichlet_l  = .FALSE.
1503                bc_dirichlet_r  = .TRUE.
1504                bc_radiation_l = .TRUE.
1505                bc_radiation_r = .FALSE.
1506             ELSEIF ( child_domain  .OR.  nesting_offline )  THEN
1507                bc_dirichlet_l = .TRUE.
1508                bc_dirichlet_r = .TRUE.
[1575]1509             ENDIF
1510
1511             IF ( bc_ns_dirrad )  THEN
[3182]1512                bc_dirichlet_n  = .TRUE.
1513                bc_dirichlet_s  = .FALSE.
1514                bc_radiation_n = .FALSE.
1515                bc_radiation_s = .TRUE.
[1575]1516             ELSEIF ( bc_ns_raddir )  THEN
[3182]1517                bc_dirichlet_n  = .FALSE.
1518                bc_dirichlet_s  = .TRUE.
1519                bc_radiation_n = .TRUE.
1520                bc_radiation_s = .FALSE.
1521             ELSEIF ( child_domain  .OR.  nesting_offline)  THEN
1522                bc_dirichlet_s = .TRUE.
1523                bc_dirichlet_n = .TRUE.
[1575]1524             ENDIF
1525
1526             DEALLOCATE( f2_sub )
1527
1528          ELSE
1529
[1931]1530             CALL restrict( f2, r )
[1575]1531
1532             ind_even_odd = even_odd_level(grid_level)  ! must be after restrict
1533
1534          ENDIF
1535
1536          p2 = 0.0_wp
1537
1538!
1539!--       Repeat the same procedure till the coarsest grid is reached
[1931]1540          CALL next_mg_level( f2, p2, p3, r )
[1575]1541
1542       ENDIF
1543
1544!
1545!--    Now follows the prolongation
1546       IF ( grid_level >= 2 )  THEN
1547
1548!
1549!--       Prolongation of the new residual. The values are transferred
1550!--       from the coarse to the next finer grid.
1551          IF ( grid_level == mg_switch_to_pe0_level+1 )  THEN
1552
1553#if defined( __parallel )
1554!
1555!--          At this level, the new residual first has to be scattered from
1556!--          PE0 to the other PEs
1557             ALLOCATE( p2_sub(nzb:mg_loc_ind(5,myid)+1,             &
1558                       mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1,   &
1559                       mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) )
1560
[1931]1561             CALL mg_scatter( p2, p2_sub )
[1575]1562
1563!
1564!--          Therefore, indices of the previous level have to be changed to
1565!--          subdomain values in between (otherwise, the prolong routine would
1566!--          expect the gathered array)
1567             nxl_mg_save = nxl_mg(grid_level-1)
1568             nxr_mg_save = nxr_mg(grid_level-1)
1569             nys_mg_save = nys_mg(grid_level-1)
1570             nyn_mg_save = nyn_mg(grid_level-1)
1571             nzt_mg_save = nzt_mg(grid_level-1)
1572             nxl_mg(grid_level-1) = mg_loc_ind(1,myid)
1573             nxr_mg(grid_level-1) = mg_loc_ind(2,myid)
1574             nys_mg(grid_level-1) = mg_loc_ind(3,myid)
1575             nyn_mg(grid_level-1) = mg_loc_ind(4,myid)
1576             nzt_mg(grid_level-1) = mg_loc_ind(5,myid)
1577
1578!
1579!--          Set switch for routine exchange_horiz, that ghostpoint exchange
1580!--          has to be carried again out from now on
1581             mg_switch_to_pe0 = .FALSE.
1582
1583!
[2021]1584!--          For non-cyclic lateral boundary conditions and in case of nesting,
1585!--          restore the in-/outflow conditions.
[3182]1586             bc_dirichlet_l = .FALSE.;  bc_dirichlet_r = .FALSE.
1587             bc_dirichlet_n = .FALSE.;  bc_dirichlet_s = .FALSE.
1588             bc_radiation_l = .FALSE.;  bc_radiation_r = .FALSE.
1589             bc_radiation_n = .FALSE.;  bc_radiation_s = .FALSE.
[1575]1590
1591             IF ( pleft == MPI_PROC_NULL )  THEN
[3182]1592                IF ( bc_lr_dirrad  .OR.  child_domain  .OR.  nesting_offline ) &
1593                THEN
1594                   bc_dirichlet_l = .TRUE.
[1575]1595                ELSEIF ( bc_lr_raddir )  THEN
[3182]1596                   bc_radiation_l = .TRUE.
[1575]1597                ENDIF
1598             ENDIF
1599
1600             IF ( pright == MPI_PROC_NULL )  THEN
1601                IF ( bc_lr_dirrad )  THEN
[3182]1602                   bc_radiation_r = .TRUE.
1603                ELSEIF ( bc_lr_raddir  .OR.  child_domain  .OR.                &
1604                         nesting_offline )  THEN
1605                   bc_dirichlet_r = .TRUE.
[1575]1606                ENDIF
1607             ENDIF
1608
1609             IF ( psouth == MPI_PROC_NULL )  THEN
1610                IF ( bc_ns_dirrad )  THEN
[3182]1611                   bc_radiation_s = .TRUE.
1612                ELSEIF ( bc_ns_raddir  .OR.  child_domain  .OR.                &
1613                         nesting_offline )  THEN
1614                   bc_dirichlet_s = .TRUE.
[1575]1615                ENDIF
1616             ENDIF
1617
1618             IF ( pnorth == MPI_PROC_NULL )  THEN
[3182]1619                IF ( bc_ns_dirrad  .OR.  child_domain  .OR.  nesting_offline ) &
1620                THEN
1621                   bc_dirichlet_n = .TRUE.
[1575]1622                ELSEIF ( bc_ns_raddir )  THEN
[3182]1623                   bc_radiation_n = .TRUE.
[1575]1624                ENDIF
1625             ENDIF
1626
[1931]1627             CALL prolong( p2_sub, p3 )
[1575]1628
1629!
1630!--          Restore the correct indices of the previous level
1631             nxl_mg(grid_level-1) = nxl_mg_save
1632             nxr_mg(grid_level-1) = nxr_mg_save
1633             nys_mg(grid_level-1) = nys_mg_save
1634             nyn_mg(grid_level-1) = nyn_mg_save
1635             nzt_mg(grid_level-1) = nzt_mg_save
1636
1637             DEALLOCATE( p2_sub )
1638#endif
1639
1640          ELSE
1641
[1931]1642             CALL prolong( p2, p3 )
[1575]1643
1644          ENDIF
1645
1646!
1647!--       Computation of the new pressure correction. Therefore,
1648!--       values from prior grids are added up automatically stage by stage.
1649          DO  i = nxl_mg(grid_level)-1, nxr_mg(grid_level)+1
1650             DO  j = nys_mg(grid_level)-1, nyn_mg(grid_level)+1
1651                DO  k = nzb, nzt_mg(grid_level)+1
1652                   p_mg(k,j,i) = p_mg(k,j,i) + p3(k,j,i)
1653                ENDDO
1654             ENDDO
1655          ENDDO
1656
1657!
1658!--       Relaxation of the new solution
[1931]1659          CALL redblack( f_mg, p_mg )
[1575]1660
1661       ENDIF
1662
1663
1664!
1665!--    The following few lines serve the steering of the multigrid scheme
1666       IF ( grid_level == maximum_grid_level )  THEN
1667
1668          GOTO 20
1669
1670       ELSEIF ( grid_level /= maximum_grid_level  .AND.  grid_level /= 1  .AND. &
1671                grid_level_count(grid_level) /= gamma_mg )  THEN
1672
1673          GOTO 10
1674
1675       ENDIF
1676
1677!
[1931]1678!--    Reset counter for the next call of poismg
[1575]1679       grid_level_count(grid_level) = 0
1680
1681!
1682!--    Continue with the next finer level. nxl..nzt have to be
1683!--    set to the finer grid values, because these variables are needed for the
1684!--    exchange of ghost points in routine exchange_horiz
1685       grid_level = grid_level + 1
1686       ind_even_odd = even_odd_level(grid_level)
1687
1688       nxl = nxl_mg(grid_level)
1689       nxr = nxr_mg(grid_level)
1690       nys = nys_mg(grid_level)
1691       nyn = nyn_mg(grid_level)
1692       nzt = nzt_mg(grid_level)
1693
1694    20 CONTINUE
1695
[1931]1696    END SUBROUTINE next_mg_level
[1575]1697
1698
1699!------------------------------------------------------------------------------!
1700! Description:
1701! ------------
[1682]1702!> Initial settings for sorting k-dimension from sequential order (alternate
1703!> even/odd) into blocks of even and odd or vice versa
[1575]1704!------------------------------------------------------------------------------!
[1682]1705    SUBROUTINE init_even_odd_blocks
[1575]1706
[1682]1707
[1575]1708       USE arrays_3d,                                                          &
[2232]1709           ONLY:  f1_mg, f2_mg, f3_mg
[1575]1710
1711       USE control_parameters,                                                 &
[1898]1712           ONLY:  grid_level, maximum_grid_level
[1575]1713
1714       USE indices,                                                            &
1715           ONLY:  nzb, nzt, nzt_mg
1716
1717       USE indices,                                                            &
[3241]1718           ONLY:  nzb, nzt_mg
[1575]1719
1720       IMPLICIT NONE
1721!
1722!--    Local variables
[1898]1723       INTEGER(iwp) ::  i     !< 
[1682]1724       INTEGER(iwp) ::  l     !<
[1575]1725
1726       LOGICAL, SAVE ::  lfirst = .TRUE.
1727
1728
1729       IF ( .NOT. lfirst )  RETURN
1730
1731       ALLOCATE( even_odd_level(maximum_grid_level) )
1732
1733       ALLOCATE( f1_mg_b(nzb:nzt+1,maximum_grid_level),                        &
1734                 f2_mg_b(nzb:nzt+1,maximum_grid_level),                        &
[2232]1735                 f3_mg_b(nzb:nzt+1,maximum_grid_level) )
[1575]1736
1737!
1738!--    Set border index between the even and odd block
1739       DO  i = maximum_grid_level, 1, -1
1740          even_odd_level(i) = nzt_mg(i) / 2
1741       ENDDO
1742
1743!
1744!--    Sort grid coefficients used in red/black scheme and for calculating the
1745!--    residual to block (even/odd) structure
1746       DO  l = maximum_grid_level, 1 , -1
1747          CALL sort_k_to_even_odd_blocks( f1_mg(nzb+1:nzt_mg(grid_level),l),   &
1748                                          f1_mg_b(nzb:nzt_mg(grid_level)+1,l), &
1749                                          l )
1750          CALL sort_k_to_even_odd_blocks( f2_mg(nzb+1:nzt_mg(grid_level),l),   &
1751                                          f2_mg_b(nzb:nzt_mg(grid_level)+1,l), &
1752                                          l )
1753          CALL sort_k_to_even_odd_blocks( f3_mg(nzb+1:nzt_mg(grid_level),l),   &
1754                                          f3_mg_b(nzb:nzt_mg(grid_level)+1,l), &
1755                                          l )
1756       ENDDO
1757
1758       lfirst = .FALSE.
1759
1760     END SUBROUTINE init_even_odd_blocks
1761
1762
1763!------------------------------------------------------------------------------!
1764! Description:
1765! ------------
[1682]1766!> Special exchange_horiz subroutine for use in redblack. Transfers only
1767!> "red" or "black" data points.
[1575]1768!------------------------------------------------------------------------------!
[3241]1769     SUBROUTINE special_exchange_horiz( p_mg, color )
[1575]1770
[1682]1771
[1575]1772       USE control_parameters,                                                 &
[4429]1773           ONLY:  grid_level
[1575]1774
[4429]1775#if defined( __parallel )
1776       USE control_parameters,                                                 &
1777           ONLY:  mg_switch_to_pe0_level, synchronous_exchange
1778#endif
1779
[1575]1780       USE indices,                                                            &
[4429]1781           ONLY:  nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
[1575]1782
[4429]1783#if defined( __parallel )
1784       USE indices,                                                            &
1785           ONLY:  nxl, nxr, nys, nyn, nzt
1786#endif
1787
[1575]1788       IMPLICIT NONE
1789
1790       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
1791                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
[1904]1792                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::       &
1793                                    p_mg   !< treated array
[1575]1794
[4429]1795       INTEGER(iwp), INTENT(IN) ::  color  !< flag for grid point type (red or black)
1796
1797#if defined ( __parallel )
[1575]1798!
1799!--    Local variables
[1904]1800       INTEGER(iwp) ::  i        !< index variable along x
1801       INTEGER(iwp) ::  i1       !< index variable along x on coarse level
1802       INTEGER(iwp) ::  i2       !< index variable along x on coarse level
[1575]1803
[1904]1804       INTEGER(iwp) ::  j        !< index variable along y
1805       INTEGER(iwp) ::  j1       !< index variable along y on coarse level
1806       INTEGER(iwp) ::  j2       !< index variable along y on coarse level
1807       INTEGER(iwp) ::  k        !< index variable along z
1808       INTEGER(iwp) ::  l        !< short for grid level
1809       INTEGER(iwp) ::  jys      !< index for lower local PE boundary along y
1810       INTEGER(iwp) ::  jyn      !< index for upper local PE boundary along y
1811       INTEGER(iwp) ::  ixl      !< index for lower local PE boundary along x
1812       INTEGER(iwp) ::  ixr      !< index for upper local PE boundary along x
[1575]1813
[1904]1814       LOGICAL      ::  synchronous_exchange_save    !< dummy to reset synchronous_exchange to prescribed value
1815
1816       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  temp  !< temporary array on next coarser grid level
1817
[1575]1818       synchronous_exchange_save   = synchronous_exchange
1819       synchronous_exchange        = .FALSE.
1820
1821       l = grid_level
1822
1823       ind_even_odd = even_odd_level(grid_level)
1824
[1904]1825!
1826!--    Restricted transfer only on finer levels with enough data.
1827!--    Restricted transfer is not possible for levels smaller or equal to
1828!--    'switch to PE0 levels', since array bounds does not fit. Moreover,
1829!--    it is not possible for the coarsest grid level, since the dimensions
1830!--    of temp are not defined. For such cases, normal exchange_horiz is called.
1831       IF ( l > 1 .AND. l > mg_switch_to_pe0_level + 1 .AND.                   &
1832            ( ngp_xz(grid_level) >= 900 .OR. ngp_yz(grid_level) >= 900 ) )  THEN
[1575]1833
[1904]1834          jys = nys_mg(grid_level-1)
1835          jyn = nyn_mg(grid_level-1)
1836          ixl = nxl_mg(grid_level-1)
1837          ixr = nxr_mg(grid_level-1)
1838          ALLOCATE( temp(nzb:nzt_mg(l-1)+1,jys-1:jyn+1,ixl-1:ixr+1) )
[1575]1839!
1840!--       Handling the even k Values
1841!--       Collecting data for the north - south exchange
1842!--       Since only every second value has to be transfered, data are stored
1843!--       on the next coarser grid level, because the arrays on that level
1844!--       have just the required size
1845          i1 = nxl_mg(grid_level-1)
1846          i2 = nxl_mg(grid_level-1)
1847
1848          DO  i = nxl_mg(l), nxr_mg(l), 2
1849             DO  j = nys_mg(l) + 2 - color, nyn_mg(l), 2
1850
1851                IF ( j == nys_mg(l) )  THEN
1852                   !DIR$ IVDEP
1853                   DO  k = ind_even_odd+1, nzt_mg(l)
1854                      temp(k-ind_even_odd,jys,i1) = p_mg(k,j,i)
1855                   ENDDO
1856                   i1 = i1 + 1
1857
1858                ENDIF
1859
1860                IF ( j == nyn_mg(l) )  THEN
1861                   !DIR$ IVDEP
1862                   DO  k = ind_even_odd+1, nzt_mg(l)
1863                      temp(k-ind_even_odd,jyn,i2) = p_mg(k,j,i)
1864                   ENDDO
1865                   i2 = i2 + 1
1866
1867                ENDIF
1868
1869             ENDDO
1870          ENDDO
1871
1872          DO  i = nxl_mg(l)+1, nxr_mg(l), 2
1873             DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
1874
1875                IF ( j == nys_mg(l) )  THEN
1876                   !DIR$ IVDEP
1877                   DO  k = ind_even_odd+1, nzt_mg(l)
1878                      temp(k-ind_even_odd,jys,i1) = p_mg(k,j,i)
1879                   ENDDO
1880                   i1 = i1 + 1
1881
1882                ENDIF
1883
1884                IF ( j == nyn_mg(l) )  THEN
1885                   !DIR$ IVDEP
1886                   DO  k = ind_even_odd+1, nzt_mg(l)
1887                      temp(k-ind_even_odd,jyn,i2) = p_mg(k,j,i)
1888                   ENDDO
1889                   i2 = i2 + 1
1890
1891                ENDIF
1892
1893             ENDDO
1894          ENDDO
1895
1896          grid_level = grid_level-1
1897
1898          nxl = nxl_mg(grid_level)
1899          nys = nys_mg(grid_level)
1900          nxr = nxr_mg(grid_level)
1901          nyn = nyn_mg(grid_level)
1902          nzt = nzt_mg(grid_level)
1903
1904          send_receive = 'ns'
1905          CALL exchange_horiz( temp, 1 )
1906
1907          grid_level = grid_level+1
1908
1909          i1 = nxl_mg(grid_level-1)
1910          i2 = nxl_mg(grid_level-1)
1911
1912          DO  i = nxl_mg(l), nxr_mg(l), 2
1913             DO  j = nys_mg(l) + 2 - color, nyn_mg(l), 2
1914
1915                IF ( j == nys_mg(l) )  THEN
1916                   !DIR$ IVDEP
1917                   DO  k = ind_even_odd+1, nzt_mg(l)
1918                      p_mg(k,nyn_mg(l)+1,i) = temp(k-ind_even_odd,jyn+1,i1)
1919                   ENDDO
1920                   i1 = i1 + 1
1921
1922                ENDIF
1923
1924                IF ( j == nyn_mg(l) )  THEN
1925                   !DIR$ IVDEP
1926                   DO  k = ind_even_odd+1, nzt_mg(l)
1927                      p_mg(k,nys_mg(l)-1,i) = temp(k-ind_even_odd,jys-1,i2)
1928                   ENDDO
1929                   i2 = i2 + 1
1930
1931                ENDIF
1932
1933             ENDDO
1934          ENDDO
1935
1936          DO  i = nxl_mg(l)+1, nxr_mg(l), 2
1937             DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
1938
1939                IF ( j == nys_mg(l) )  THEN
1940                   !DIR$ IVDEP
1941                   DO  k = ind_even_odd+1, nzt_mg(l)
1942                      p_mg(k,nyn_mg(l)+1,i) = temp(k-ind_even_odd,jyn+1,i1)
1943                   ENDDO
1944                   i1 = i1 + 1
1945
1946                ENDIF
1947
1948                IF ( j == nyn_mg(l) )  THEN
1949                   !DIR$ IVDEP
1950                   DO  k = ind_even_odd+1, nzt_mg(l)
1951                      p_mg(k,nys_mg(l)-1,i) = temp(k-ind_even_odd,jys-1,i2)
1952                   ENDDO
1953                   i2 = i2 + 1
1954
1955                ENDIF
1956
1957             ENDDO
1958          ENDDO
1959
1960!
1961!--       Collecting data for the left - right exchange
1962!--       Since only every second value has to be transfered, data are stored
1963!--       on the next coarser grid level, because the arrays on that level
1964!--       have just the required size
1965          j1 = nys_mg(grid_level-1)
1966          j2 = nys_mg(grid_level-1)
1967
1968          DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
1969             DO  i = nxl_mg(l), nxr_mg(l), 2
1970
1971                IF ( i == nxl_mg(l) )  THEN
1972                   !DIR$ IVDEP
1973                   DO  k = ind_even_odd+1, nzt_mg(l)
1974                      temp(k-ind_even_odd,j1,ixl) = p_mg(k,j,i)
1975                   ENDDO
1976                   j1 = j1 + 1
1977
1978                ENDIF
1979
1980                IF ( i == nxr_mg(l) )  THEN
1981                   !DIR$ IVDEP
1982                   DO  k = ind_even_odd+1, nzt_mg(l)
1983                      temp(k-ind_even_odd,j2,ixr) = p_mg(k,j,i)
1984                   ENDDO
1985                   j2 = j2 + 1
1986
1987                ENDIF
1988
1989             ENDDO
1990          ENDDO
1991
1992          DO j = nys_mg(l) + (color-1), nyn_mg(l), 2
1993             DO  i = nxl_mg(l)+1, nxr_mg(l), 2
1994
1995                IF ( i == nxl_mg(l) )  THEN
1996                   !DIR$ IVDEP
1997                   DO  k = ind_even_odd+1, nzt_mg(l)
1998                      temp(k-ind_even_odd,j1,ixl) = p_mg(k,j,i)
1999                   ENDDO
2000                   j1 = j1 + 1
2001
2002                ENDIF
2003
2004                IF ( i == nxr_mg(l) )  THEN
2005                   !DIR$ IVDEP
2006                   DO  k = ind_even_odd+1, nzt_mg(l)
2007                      temp(k-ind_even_odd,j2,ixr) = p_mg(k,j,i)
2008                   ENDDO
2009                   j2 = j2 + 1
2010
2011                ENDIF
2012
2013             ENDDO
2014          ENDDO
2015
2016          grid_level = grid_level-1
2017          send_receive = 'lr'
2018
2019          CALL exchange_horiz( temp, 1 )
2020
2021          grid_level = grid_level+1
2022
2023          j1 = nys_mg(grid_level-1)
2024          j2 = nys_mg(grid_level-1)
2025
2026          DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
2027             DO  i = nxl_mg(l), nxr_mg(l), 2
2028
2029                IF ( i == nxl_mg(l) )  THEN
2030                   !DIR$ IVDEP
2031                   DO  k = ind_even_odd+1, nzt_mg(l)
2032                      p_mg(k,j,nxr_mg(l)+1)  = temp(k-ind_even_odd,j1,ixr+1)
2033                   ENDDO
2034                   j1 = j1 + 1
2035
2036                ENDIF
2037
2038                IF ( i == nxr_mg(l) )  THEN
2039                   !DIR$ IVDEP
2040                   DO  k = ind_even_odd+1, nzt_mg(l)
2041                      p_mg(k,j,nxl_mg(l)-1)  = temp(k-ind_even_odd,j2,ixl-1)
2042                   ENDDO
2043                   j2 = j2 + 1
2044
2045                ENDIF
2046
2047             ENDDO
2048          ENDDO
2049
2050          DO j = nys_mg(l) + (color-1), nyn_mg(l), 2
2051             DO  i = nxl_mg(l)+1, nxr_mg(l), 2
2052
2053                IF ( i == nxl_mg(l) )  THEN
2054                   !DIR$ IVDEP
2055                   DO  k = ind_even_odd+1, nzt_mg(l)
2056                      p_mg(k,j,nxr_mg(l)+1)  = temp(k-ind_even_odd,j1,ixr+1)
2057                   ENDDO
2058                   j1 = j1 + 1
2059
2060                ENDIF
2061
2062                IF ( i == nxr_mg(l) )  THEN
2063                   !DIR$ IVDEP
2064                   DO  k = ind_even_odd+1, nzt_mg(l)
2065                      p_mg(k,j,nxl_mg(l)-1)  = temp(k-ind_even_odd,j2,ixl-1)
2066                   ENDDO
2067                   j2 = j2 + 1
2068
2069                ENDIF
2070
2071             ENDDO
2072          ENDDO
2073
2074!
2075!--       Now handling the even k values
2076!--       Collecting data for the north - south exchange
2077!--       Since only every second value has to be transfered, data are stored
2078!--       on the next coarser grid level, because the arrays on that level
2079!--       have just the required size
2080          i1 = nxl_mg(grid_level-1)
2081          i2 = nxl_mg(grid_level-1)
2082
2083          DO  i = nxl_mg(l), nxr_mg(l), 2
2084             DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
2085
2086                IF ( j == nys_mg(l) )  THEN
2087                   !DIR$ IVDEP
2088                   DO  k = nzb+1, ind_even_odd
2089                      temp(k,jys,i1) = p_mg(k,j,i)
2090                   ENDDO
2091                   i1 = i1 + 1
2092
2093                ENDIF
2094
2095                IF ( j == nyn_mg(l) )  THEN
2096                   !DIR$ IVDEP
2097                   DO  k = nzb+1, ind_even_odd
2098                      temp(k,jyn,i2) = p_mg(k,j,i)
2099                   ENDDO
2100                   i2 = i2 + 1
2101
2102                ENDIF
2103
2104             ENDDO
2105          ENDDO
2106
2107          DO  i = nxl_mg(l)+1, nxr_mg(l), 2
2108             DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
2109
2110                IF ( j == nys_mg(l) )  THEN
2111                   !DIR$ IVDEP
2112                   DO  k = nzb+1, ind_even_odd
2113                      temp(k,jys,i1) = p_mg(k,j,i)
2114                   ENDDO
2115                   i1 = i1 + 1
2116
2117                ENDIF
2118
2119                IF ( j == nyn_mg(l) )  THEN
2120                   !DIR$ IVDEP
2121                   DO  k = nzb+1, ind_even_odd
2122                      temp(k,jyn,i2) = p_mg(k,j,i)
2123                   ENDDO
2124                   i2 = i2 + 1
2125
2126                ENDIF
2127
2128             ENDDO
2129          ENDDO
2130
2131          grid_level = grid_level-1
2132
2133          send_receive = 'ns'
2134          CALL exchange_horiz( temp, 1 )
2135
2136          grid_level = grid_level+1
2137
2138          i1 = nxl_mg(grid_level-1)
2139          i2 = nxl_mg(grid_level-1)
2140
2141          DO  i = nxl_mg(l), nxr_mg(l), 2
2142             DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
2143
2144                IF ( j == nys_mg(l) )  THEN
2145                   !DIR$ IVDEP
2146                   DO  k = nzb+1, ind_even_odd
2147                      p_mg(k,nyn_mg(l)+1,i) = temp(k,jyn+1,i1)
2148                   ENDDO
2149                   i1 = i1 + 1
2150
2151                ENDIF
2152
2153                IF ( j == nyn_mg(l) )  THEN
2154                   !DIR$ IVDEP
2155                   DO  k = nzb+1, ind_even_odd
2156                      p_mg(k,nys_mg(l)-1,i) = temp(k,jys-1,i2)
2157                   ENDDO
2158                   i2 = i2 + 1
2159
2160                ENDIF
2161
2162             ENDDO
2163          ENDDO
2164
2165          DO  i = nxl_mg(l)+1, nxr_mg(l), 2
2166             DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
2167
2168                IF ( j == nys_mg(l) )  THEN
2169                   !DIR$ IVDEP
2170                   DO  k = nzb+1, ind_even_odd
2171                      p_mg(k,nyn_mg(l)+1,i) = temp(k,jyn+1,i1)
2172                   ENDDO
2173                   i1 = i1 + 1
2174
2175                ENDIF
2176
2177                IF ( j == nyn_mg(l) )  THEN
2178                   !DIR$ IVDEP
2179                   DO  k = nzb+1, ind_even_odd
2180                      p_mg(k,nys_mg(l)-1,i) = temp(k,jys-1,i2)
2181                   ENDDO
2182                   i2 = i2 + 1
2183
2184                ENDIF
2185
2186             ENDDO
2187          ENDDO
2188
2189          j1 = nys_mg(grid_level-1)
2190          j2 = nys_mg(grid_level-1)
2191
2192          DO  i = nxl_mg(l), nxr_mg(l), 2
2193             DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
2194
2195                IF ( i == nxl_mg(l) )  THEN
2196                   !DIR$ IVDEP
2197                   DO  k = nzb+1, ind_even_odd
2198                      temp(k,j1,ixl) = p_mg(k,j,i)
2199                   ENDDO
2200                   j1 = j1 + 1
2201
2202                ENDIF
2203
2204                IF ( i == nxr_mg(l) )  THEN
2205                   !DIR$ IVDEP
2206                   DO  k = nzb+1, ind_even_odd
2207                      temp(k,j2,ixr) = p_mg(k,j,i)
2208                   ENDDO
2209                   j2 = j2 + 1
2210
2211                ENDIF
2212
2213             ENDDO
2214          ENDDO
2215
2216          DO  i = nxl_mg(l)+1, nxr_mg(l), 2
2217             DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
2218
2219                IF ( i == nxl_mg(l) )  THEN
2220                   !DIR$ IVDEP
2221                   DO  k = nzb+1, ind_even_odd
2222                      temp(k,j1,ixl) = p_mg(k,j,i)
2223                   ENDDO
2224                   j1 = j1 + 1
2225
2226                ENDIF
2227
2228                IF ( i == nxr_mg(l) )  THEN
2229                   !DIR$ IVDEP
2230                   DO  k = nzb+1, ind_even_odd
2231                      temp(k,j2,ixr) = p_mg(k,j,i)
2232                   ENDDO
2233                   j2 = j2 + 1
2234
2235                ENDIF
2236
2237             ENDDO
2238          ENDDO
2239
2240          grid_level = grid_level-1
2241
2242          send_receive = 'lr'
2243          CALL exchange_horiz( temp, 1 )
2244
2245          grid_level = grid_level+1
2246
2247          nxl = nxl_mg(grid_level)
2248          nys = nys_mg(grid_level)
2249          nxr = nxr_mg(grid_level)
2250          nyn = nyn_mg(grid_level)
2251          nzt = nzt_mg(grid_level)
2252
2253          j1 = nys_mg(grid_level-1)
2254          j2 = nys_mg(grid_level-1)
2255
2256          DO  i = nxl_mg(l), nxr_mg(l), 2
2257             DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
2258
2259                IF ( i == nxl_mg(l) )  THEN
2260                   !DIR$ IVDEP
2261                   DO  k = nzb+1, ind_even_odd
2262                      p_mg(k,j,nxr_mg(l)+1)  = temp(k,j1,ixr+1)
2263                   ENDDO
2264                   j1 = j1 + 1
2265
2266                ENDIF
2267
2268                IF ( i == nxr_mg(l) )  THEN
2269                   !DIR$ IVDEP
2270                   DO  k = nzb+1, ind_even_odd
2271                      p_mg(k,j,nxl_mg(l)-1)  = temp(k,j2,ixl-1)
2272                   ENDDO
2273                   j2 = j2 + 1
2274
2275                ENDIF
2276
2277             ENDDO
2278          ENDDO
2279
2280          DO  i = nxl_mg(l)+1, nxr_mg(l), 2
2281             DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
2282
2283                IF ( i == nxl_mg(l) )  THEN
[4432]2284                   !DIR$ IVDEP
[1575]2285                   DO  k = nzb+1, ind_even_odd
2286                      p_mg(k,j,nxr_mg(l)+1)  = temp(k,j1,ixr+1)
2287                   ENDDO
2288                   j1 = j1 + 1
2289
2290                ENDIF
2291
2292                IF ( i == nxr_mg(l) )  THEN
2293                   !DIR$ IVDEP
2294                   DO  k = nzb+1, ind_even_odd
2295                      p_mg(k,j,nxl_mg(l)-1)  = temp(k,j2,ixl-1)
2296                   ENDDO
2297                   j2 = j2 + 1
2298
2299                ENDIF
2300
2301             ENDDO
2302          ENDDO
2303
[1904]2304          DEALLOCATE( temp )
2305
[1575]2306       ELSE
[1609]2307
[1575]2308!
2309!--       Standard horizontal ghost boundary exchange for small coarse grid
2310!--       levels, where the transfer time is latency bound
2311          CALL exchange_horiz( p_mg, 1 )
2312
2313       ENDIF
2314
2315!
2316!--    Reset values to default PALM setup
2317       synchronous_exchange   = synchronous_exchange_save
2318       send_receive = 'al'
[1609]2319#else
[1575]2320
[1609]2321!
[4429]2322!--    Next line is to avoid compile error due to unused dummy argument
2323       IF ( color == 1234567 )  RETURN
2324!
[1609]2325!--    Standard horizontal ghost boundary exchange for small coarse grid
2326!--    levels, where the transfer time is latency bound
2327       CALL exchange_horiz( p_mg, 1 )
2328#endif
2329
[1575]2330    END SUBROUTINE special_exchange_horiz
2331
2332 END MODULE poismg_mod
Note: See TracBrowser for help on using the repository browser.