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

Last change on this file since 3689 was 3655, checked in by knoop, 6 years ago

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

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