source: palm/trunk/SOURCE/init_pegrid.f90 @ 4069

Last change on this file since 4069 was 4045, checked in by raasch, 5 years ago

bugfix: kind attribute added to NINT function to allow for large integers which may appear in case of default recycling width and small grid spacings

  • Property svn:keywords set to Id
File size: 56.0 KB
Line 
1!> @file init_pegrid.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2019 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: init_pegrid.f90 4045 2019-06-21 10:58:47Z Giersch $
27! bugfix: kind attribute added to nint function to allow for large integers which may appear in
28! case of default recycling width and small grid spacings
29!
30! 3999 2019-05-23 16:09:37Z suehring
31! Spend 3 ghost points also in case of pw-scheme when nesting is applied
32!
33! 3897 2019-04-15 11:51:14Z suehring
34! Minor revision of multigrid check; give warning instead of an abort.
35!
36! 3890 2019-04-12 15:59:20Z suehring
37! Check if grid coarsening is possible on subdomain, in order to avoid that
38! multigrid approach effectively reduces to a Gauss-Seidel scheme.
39!
40! 3885 2019-04-11 11:29:34Z kanani
41! Changes related to global restructuring of location messages and introduction
42! of additional debug messages
43!
44! 3884 2019-04-10 13:31:55Z Giersch
45! id_recycling is only calculated in case of tubulent inflow
46!
47! 3761 2019-02-25 15:31:42Z raasch
48! unused variable removed
49!
50! 3655 2019-01-07 16:51:22Z knoop
51! variables documented
52!
53! 3552 2018-11-22 10:28:35Z suehring
54! Introduce new MPI-datatype for ghostpoint exchange of 2D 8-bit Integer arrays
55!
56! 3542 2018-11-20 17:04:13Z suehring
57! Bugfix in setting number of ghost layers in neutral case
58!
59! 3341 2018-10-15 10:31:27Z suehring
60! unused variables removed
61!
62! 3183 2018-07-27 14:25:55Z suehring
63! Rename variables and boundary conditions in mesoscale-offline nesting mode
64!
65! 3182 2018-07-27 13:36:03Z suehring
66! bugfix: wrong error number in r3057 revised
67!
68! 3057 2018-06-05 09:03:41Z raasch
69! bugfix: check that nz is even in case that optimized multigrid is used
70!
71! 3049 2018-05-29 13:52:36Z Giersch
72! Error messages revised
73!
74! 3045 2018-05-28 07:55:41Z Giersch
75! Error messages revised
76!
77! 2938 2018-03-27 15:52:42Z suehring
78! - No checks for domain decomposition in case of turbulence generator
79!  (is done in stg module)
80! - Introduce ids to indicate lateral processors for turbulence generator
81!
82! 2936 2018-03-27 14:49:27Z suehring
83! Variable use_synthetic_turbulence_generator has been abbreviated
84!
85! 2718 2018-01-02 08:49:38Z maronga
86! Corrected "Former revisions" section
87!
88! 2696 2017-12-14 17:12:51Z kanani
89! Change in file header (GPL part)
90! 3D-Integer exchange on multigrid level (MS)
91! Forcing implemented (MS)
92!
93! 2600 2017-11-01 14:11:20Z raasch
94! calculation of block-I/O quantitites removed (is now done in parin)
95!
96! 2516 2017-10-04 11:03:04Z suehring
97! Remove tabs
98!
99! 2514 2017-10-04 09:52:37Z suehring
100! Redundant preprocessor directives removed
101!
102! 2372 2017-08-25 12:37:32Z sward
103! Shifted cyclic boundary conditions implemented
104!
105! 2365 2017-08-21 14:59:59Z kanani
106! Vertical nesting implemented (SadiqHuq)
107!
108! 2300 2017-06-29 13:31:14Z raasch
109! host-specific settings removed
110!
111! 2298 2017-06-29 09:28:18Z raasch
112! MPI2 related parts removed
113!
114! 2271 2017-06-09 12:34:55Z sward
115! Error message changed
116!
117! 2259 2017-06-08 09:09:11Z gronemeier
118! Implemented synthetic turbulence generator
119!
120! 2238 2017-05-31 16:49:16Z suehring
121! Remove unnecessary module load of pmc_interface
122!
123! 2231 2017-05-30 16:44:33Z suehring
124!
125! 2200 2017-04-11 11:37:51Z suehring
126! monotonic_adjustment removed
127!
128! 2197 2017-03-24 02:25:00Z raasch
129! bugfix: do not allow odd values for nz at the coarsest grid level in case of
130! optimized multigrid solver
131!
132! 2180 2017-03-17 13:33:05Z hellstea
133! Checks to ensure (2178) that pdims match the grid dimensions in the
134! automatic determination of pdims are canceled as unnecessary
135!
136! 2178 2017-03-17 11:07:39Z hellstea
137! Checks to ensure that pdims match the grid dimensions are added in the
138! automatic determination of pdims
139!
140! 2050 2016-11-08 15:00:55Z gronemeier
141! Implement turbulent outflow condition
142!
143! 2000 2016-08-20 18:09:15Z knoop
144! Forced header and separation lines into 80 columns
145!
146! 1968 2016-07-18 12:01:49Z suehring
147! Extent MPI-datatypes for exchange of 2D-INTEGER arrays on coarser multigrid
148! level 
149!
150! 1964 2016-07-14 15:35:18Z hellstea
151! Bugfix: erroneous setting of nest_bound_l/r/s/n = .TRUE. for vertical nesting mode removed.
152!
153! 1923 2016-05-31 16:37:07Z boeske
154! Initial version of purely vertical nesting introduced.
155!
156! 1922 2016-05-31 16:36:08Z boeske
157! Bugfix: array transposition checks restricted to cases if a fourier
158! transform is used , removed unused variable nnx_z
159!
160! 1833 2016-04-07 14:23:03Z raasch
161! spectra related variables moved to spectra_mod
162!
163! 1815 2016-04-06 13:49:59Z raasch
164! cpp-directives for intel openmp bug removed
165!
166! 1804 2016-04-05 16:30:18Z maronga
167! Removed code for parameter file check (__check)
168!
169! 1779 2016-03-03 08:01:28Z raasch
170! changes regarding nested domain removed: virtual PE grid will be automatically
171! calculated for nested runs too
172!
173! 1764 2016-02-28 12:45:19Z raasch
174! cpp-statements for nesting removed
175!
176! 1762 2016-02-25 12:31:13Z hellstea
177! Introduction of nested domain feature
178!
179! 1682 2015-10-07 23:56:08Z knoop
180! Code annotations made doxygen readable
181!
182! 1677 2015-10-02 13:25:23Z boeske
183! New MPI-data types for exchange of 3D integer arrays.
184!
185! 1575 2015-03-27 09:56:27Z raasch
186! adjustments for psolver-queries, calculation of ngp_xz added
187!
188! 1565 2015-03-09 20:59:31Z suehring
189! Refine if-clause for setting nbgp.
190!
191! 1557 2015-03-05 16:43:04Z suehring
192! Adjustment for monotonic limiter
193!
194! 1468 2014-09-24 14:06:57Z maronga
195! Adapted for use on up to 6-digit processor cores
196!
197! 1435 2014-07-21 10:37:02Z keck
198! bugfix: added missing parameter coupling_mode_remote to ONLY-attribute
199!
200! 1402 2014-05-09 14:25:13Z raasch
201! location messages modified
202!
203! 1384 2014-05-02 14:31:06Z raasch
204! location messages added
205!
206! 1353 2014-04-08 15:21:23Z heinze
207! REAL constants provided with KIND-attribute
208!
209! 1322 2014-03-20 16:38:49Z raasch
210! REAL functions provided with KIND-attribute
211!
212! 1320 2014-03-20 08:40:49Z raasch
213! ONLY-attribute added to USE-statements,
214! kind-parameters added to all INTEGER and REAL declaration statements,
215! kinds are defined in new module kinds,
216! revision history before 2012 removed,
217! comment fields (!:) to be used for variable explanations added to
218! all variable declaration statements
219!
220! 1304 2014-03-12 10:29:42Z raasch
221! bugfix: single core MPI runs missed some settings of transpose indices
222!
223! 1212 2013-08-15 08:46:27Z raasch
224! error message for poisfft_hybrid removed
225!
226! 1159 2013-05-21 11:58:22Z fricke
227! dirichlet/neumann and neumann/dirichlet removed
228!
229! 1139 2013-04-18 07:25:03Z raasch
230! bugfix for calculating the id of the PE carrying the recycling plane
231!
232! 1111 2013-03-08 23:54:10Z raasch
233! initialization of poisfft moved to module poisfft
234!
235! 1092 2013-02-02 11:24:22Z raasch
236! unused variables removed
237!
238! 1056 2012-11-16 15:28:04Z raasch
239! Indices for arrays n.._mg start from zero due to definition of arrays f2 and
240! p2 as automatic arrays in recursive subroutine next_mg_level
241!
242! 1041 2012-11-06 02:36:29Z raasch
243! a 2d virtual processor topology is used by default for all machines
244!
245! 1036 2012-10-22 13:43:42Z raasch
246! code put under GPL (PALM 3.9)
247!
248! 1003 2012-09-14 14:35:53Z raasch
249! subdomains must have identical size (grid matching = "match" removed)
250!
251! 1001 2012-09-13 14:08:46Z raasch
252! all actions concerning upstream-spline-method removed
253!
254! 978 2012-08-09 08:28:32Z fricke
255! dirichlet/neumann and neumann/dirichlet added
256! nxlu and nysv are also calculated for inflow boundary
257!
258! 809 2012-01-30 13:32:58Z maronga
259! Bugfix: replaced .AND. and .NOT. with && and ! in the preprocessor directives
260!
261! 807 2012-01-25 11:53:51Z maronga
262! New cpp directive "__check" implemented which is used by check_namelist_files
263!
264! Revision 1.1  1997/07/24 11:15:09  raasch
265! Initial revision
266!
267!
268! Description:
269! ------------
270!> Determination of the virtual processor topology (if not prescribed by the
271!> user)and computation of the grid point number and array bounds of the local
272!> domains.
273!> @todo: remove MPI-data types for 2D exchange on coarse multigrid level (not
274!>        used any more)
275!------------------------------------------------------------------------------!
276 SUBROUTINE init_pegrid
277 
278
279    USE control_parameters,                                                    &
280        ONLY:  bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, bc_dirichlet_s, &
281               bc_lr, bc_ns, bc_radiation_l, bc_radiation_n, bc_radiation_r,   &
282               bc_radiation_s, coupling_mode, coupling_topology, gathered_size,&
283               grid_level, grid_level_count, maximum_grid_level,               &
284               message_string, mg_switch_to_pe0_level, momentum_advec,         &
285               psolver, outflow_source_plane, recycling_width, scalar_advec,   &
286               subdomain_size, turbulent_inflow, turbulent_outflow, y_shift
287
288    USE grid_variables,                                                        &
289        ONLY:  dx
290       
291    USE indices,                                                               &
292        ONLY:  mg_loc_ind, nbgp, nnx, nny, nnz, nx, nx_a, nx_o, nxl, nxl_mg,   &
293               nxlu, nxr, nxr_mg, ny, ny_a, ny_o, nyn, nyn_mg, nys, nys_mg,    &
294               nysv, nz, nzb, nzt, nzt_mg, wall_flags_1, wall_flags_2,         &
295               wall_flags_3, wall_flags_4, wall_flags_5, wall_flags_6,         &
296               wall_flags_7, wall_flags_8, wall_flags_9, wall_flags_10
297
298    USE kinds
299     
300    USE pegrid
301   
302    USE pmc_interface,                                                         &
303        ONLY:  nested_run
304     
305    USE spectra_mod,                                                           &
306        ONLY:  calculate_spectra, dt_dosp
307
308    USE synthetic_turbulence_generator_mod,                                    &
309        ONLY:  id_stg_left, id_stg_north, id_stg_right, id_stg_south,          &
310               use_syn_turb_gen
311
312    USE transpose_indices,                                                     &
313        ONLY:  nxl_y, nxl_yd, nxl_z, nxr_y, nxr_yd, nxr_z, nyn_x, nyn_z, nys_x,&
314               nys_z, nzb_x, nzb_y, nzb_yd, nzt_x, nzt_yd, nzt_y
315
316    USE vertical_nesting_mod,                                                  &
317        ONLY:  vnested, vnest_init_pegrid_domain, vnest_init_pegrid_rank
318
319    IMPLICIT NONE
320
321    INTEGER(iwp) ::  i                        !< running index over number of processors or number of multigrid level
322    INTEGER(iwp) ::  id_inflow_l              !< ID indicating processors located at the left inflow boundary
323    INTEGER(iwp) ::  id_outflow_l             !< local value of id_outflow
324    INTEGER(iwp) ::  id_outflow_source_l      !< local value of id_outflow_source
325    INTEGER(iwp) ::  id_recycling_l           !< ID indicating processors located at the recycling plane
326    INTEGER(iwp) ::  id_stg_left_l            !< left lateral boundary local core id in case of turbulence generator 
327    INTEGER(iwp) ::  id_stg_north_l           !< north lateral boundary local core id in case of turbulence generator 
328    INTEGER(iwp) ::  id_stg_right_l           !< right lateral boundary local core id in case of turbulence generator 
329    INTEGER(iwp) ::  id_stg_south_l           !< south lateral boundary local core id in case of turbulence generator 
330    INTEGER(iwp) ::  ind(5)                   !< array containing the subdomain bounds
331    INTEGER(iwp) ::  j                        !< running index, used for various loops
332    INTEGER(iwp) ::  k                        !< number of vertical grid points in different multigrid level
333    INTEGER(iwp) ::  maximum_grid_level_l     !< maximum number of grid level without switching to PE 0
334    INTEGER(iwp) ::  mg_levels_x              !< maximum number of grid level allowed along x-direction
335    INTEGER(iwp) ::  mg_levels_y              !< maximum number of grid level allowed along y-direction
336    INTEGER(iwp) ::  mg_levels_z              !< maximum number of grid level allowed along z-direction
337    INTEGER(iwp) ::  mg_switch_to_pe0_level_l !< maximum number of grid level with switching to PE 0
338    INTEGER(iwp) ::  nnx_y                    !< quotient of number of grid points along x-direction and number of PEs used along y-direction
339    INTEGER(iwp) ::  nny_x                    !< quotient of number of grid points along y-direction and number of PEs used along x-direction
340    INTEGER(iwp) ::  nny_z                    !< quotient of number of grid points along y-direction and number of PEs used along x-direction
341    INTEGER(iwp) ::  nnz_x                    !< quotient of number of grid points along z-direction and number of PEs used along x-direction
342    INTEGER(iwp) ::  nnz_y                    !< quotient of number of grid points along z-direction and number of PEs used along x-direction
343    INTEGER(iwp) ::  numproc_sqr              !< square root of the number of processors
344    INTEGER(iwp) ::  nxl_l                    !< lower index bound along x-direction on subdomain and different multigrid level
345    INTEGER(iwp) ::  nxr_l                    !< upper index bound along x-direction on subdomain and different multigrid level
346    INTEGER(iwp) ::  nyn_l                    !< lower index bound along y-direction on subdomain and different multigrid level
347    INTEGER(iwp) ::  nys_l                    !< upper index bound along y-direction on subdomain and different multigrid level
348    INTEGER(iwp) ::  nzb_l                    !< lower index bound along z-direction on subdomain and different multigrid level
349    INTEGER(iwp) ::  nzt_l                    !< upper index bound along z-direction on subdomain and different multigrid level
350!$  INTEGER(iwp) ::  omp_get_num_threads      !< number of OpenMP threads
351
352    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  ind_all !< dummy array containing index bounds on subdomain, used for gathering
353    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nxlf    !< lower index bound allong x-direction for every PE
354    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nxrf    !< upper index bound allong x-direction for every PE
355    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nynf    !< lower index bound allong y-direction for every PE
356    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nysf    !< lower index bound allong y-direction for every PE
357
358    INTEGER(iwp), DIMENSION(2) ::  pdims_remote         !< number of PEs used for coupled model (only in atmospher-ocean coupling)
359    INTEGER(iwp)               ::  lcoord(2)            !< PE coordinates of left neighbor along x and y
360    INTEGER(iwp)               ::  rcoord(2)            !< PE coordinates of right neighbor along x and y
361
362!
363!-- Get the number of OpenMP threads
364    !$OMP PARALLEL
365!$  threads_per_task = omp_get_num_threads()
366    !$OMP END PARALLEL
367
368
369#if defined( __parallel )
370
371    CALL location_message( 'creating virtual PE grids + MPI derived data types', 'start' )
372
373!
374!-- Determine the processor topology or check it, if prescribed by the user
375    IF ( npex == -1  .AND.  npey == -1 )  THEN
376
377!
378!--    Automatic determination of the topology
379       numproc_sqr = SQRT( REAL( numprocs, KIND=wp ) )
380       pdims(1)    = MAX( numproc_sqr , 1 )
381       DO  WHILE ( MOD( numprocs , pdims(1) ) /= 0 )
382          pdims(1) = pdims(1) - 1
383       ENDDO
384       pdims(2) = numprocs / pdims(1)
385
386    ELSEIF ( npex /= -1  .AND.  npey /= -1 )  THEN
387
388!
389!--    Prescribed by user. Number of processors on the prescribed topology
390!--    must be equal to the number of PEs available to the job
391       IF ( ( npex * npey ) /= numprocs )  THEN
392          WRITE( message_string, * ) 'number of PEs of the prescribed ',       &
393              'topology (', npex*npey,') does not match & the number of ',     &
394              'PEs available to the job (', numprocs, ')'
395          CALL message( 'init_pegrid', 'PA0221', 1, 2, 0, 6, 0 )
396       ENDIF
397       pdims(1) = npex
398       pdims(2) = npey
399
400    ELSE
401!
402!--    If the processor topology is prescribed by the user, the number of
403!--    PEs must be given in both directions
404       message_string = 'if the processor topology is prescribed by th' //     &
405                'e user & both values of "npex" and "npey" must be given' //   &
406                ' in the &NAMELIST-parameter file'
407       CALL message( 'init_pegrid', 'PA0222', 1, 2, 0, 6, 0 )
408
409    ENDIF
410
411!
412!-- If necessary, set horizontal boundary conditions to non-cyclic
413    IF ( bc_lr /= 'cyclic' )  cyclic(1) = .FALSE.
414    IF ( bc_ns /= 'cyclic' )  cyclic(2) = .FALSE.
415
416
417!
418!-- Create the virtual processor grid
419    CALL MPI_CART_CREATE( comm_palm, ndim, pdims, cyclic, reorder, &
420                          comm2d, ierr )
421    CALL MPI_COMM_RANK( comm2d, myid, ierr )
422    WRITE (myid_char,'(''_'',I6.6)')  myid
423
424    CALL MPI_CART_COORDS( comm2d, myid, ndim, pcoord, ierr )
425    CALL MPI_CART_SHIFT( comm2d, 0, 1, pleft, pright, ierr )
426    CALL MPI_CART_SHIFT( comm2d, 1, 1, psouth, pnorth, ierr )
427!
428!-- In case of cyclic boundary conditions, a y-shift at the boundaries in
429!-- x-direction can be introduced via parameter y_shift. The shift is done
430!-- by modifying the processor grid in such a way that processors located
431!-- at the x-boundary communicate across it to processors with y-coordinate
432!-- shifted by y_shift relative to their own. This feature can not be used
433!-- in combination with an fft pressure solver. It has been implemented to
434!-- counter the effect of streak structures in case of cyclic boundary
435!-- conditions. For a description of these see Munters
436!-- (2016; dx.doi.org/10.1063/1.4941912)
437!--
438!-- Get coordinates of left and right neighbor on PE grid
439    IF ( y_shift /= 0 ) THEN
440
441       IF ( bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic' )  THEN
442          message_string = 'y_shift /= 0 is only allowed for cyclic ' //       &
443                           'boundary conditions in both directions '
444          CALL message( 'check_parameters', 'PA0467', 1, 2, 0, 6, 0 )
445       ENDIF
446       IF ( TRIM( psolver ) /= 'multigrid' .AND.                               &
447            TRIM( psolver ) /= 'multigrid_noopt')                              &
448       THEN
449          message_string = 'y_shift /= 0 requires a multigrid pressure solver '
450          CALL message( 'check_parameters', 'PA0468', 1, 2, 0, 6, 0 )
451       ENDIF
452
453       CALL MPI_CART_COORDS( comm2d, pright, ndim, rcoord, ierr )
454       CALL MPI_CART_COORDS( comm2d, pleft, ndim, lcoord, ierr )
455
456!
457!--    If the x(y)-coordinate of the right (left) neighbor is smaller (greater)
458!--    than that of the calling process, then the calling process is located on
459!--    the right (left) boundary of the processor grid. In that case,
460!--    the y-coordinate of that neighbor is increased (decreased) by y_shift.
461!--    The rank of the process with that coordinate is then inquired and the
462!--    neighbor rank for MPI_SENDRECV, pright (pleft) is set to it.
463!--    In this way, the calling process receives a new right (left) neighbor
464!--    for all future MPI_SENDRECV calls. That neighbor has a y-coordinate
465!--    of y+(-)y_shift, where y is the original right (left) neighbor's
466!--    y-coordinate. The modulo-operation ensures that if the neighbor's
467!--    y-coordinate exceeds the grid-boundary, it will be relocated to
468!--    the opposite part of the grid cyclicly.
469       IF ( rcoord(1) < pcoord(1) ) THEN
470          rcoord(2) = MODULO( rcoord(2) + y_shift, pdims(2) )
471          CALL MPI_CART_RANK( comm2d, rcoord, pright, ierr )
472       ENDIF
473
474       IF ( lcoord(1) > pcoord(1) ) THEN
475          lcoord(2) = MODULO( lcoord(2) - y_shift, pdims(2) )
476          CALL MPI_CART_RANK( comm2d, lcoord, pleft, ierr )
477       ENDIF
478    ENDIF
479!
480!-- Vertical nesting: store four lists that identify partner ranks to exchange
481!-- data
482    IF ( vnested )  CALL vnest_init_pegrid_rank
483
484!
485!-- Determine sub-topologies for transpositions
486!-- Transposition from z to x:
487    remain_dims(1) = .TRUE.
488    remain_dims(2) = .FALSE.
489    CALL MPI_CART_SUB( comm2d, remain_dims, comm1dx, ierr )
490    CALL MPI_COMM_RANK( comm1dx, myidx, ierr )
491!
492!-- Transposition from x to y
493    remain_dims(1) = .FALSE.
494    remain_dims(2) = .TRUE.
495    CALL MPI_CART_SUB( comm2d, remain_dims, comm1dy, ierr )
496    CALL MPI_COMM_RANK( comm1dy, myidy, ierr )
497
498
499!
500!-- Calculate array bounds along x-direction for every PE.
501    ALLOCATE( nxlf(0:pdims(1)-1), nxrf(0:pdims(1)-1), nynf(0:pdims(2)-1),      &
502              nysf(0:pdims(2)-1) )
503
504    IF ( MOD( nx+1 , pdims(1) ) /= 0 )  THEN
505       WRITE( message_string, * ) 'x-direction: gridpoint number (',nx+1,') ', &
506                               'is not an& integral divisor of the number ',    &
507                               'of processors (', pdims(1),')'
508       CALL message( 'init_pegrid', 'PA0225', 1, 2, 0, 6, 0 )
509    ELSE
510       nnx  = ( nx + 1 ) / pdims(1)
511    ENDIF   
512
513!
514!-- Left and right array bounds, number of gridpoints
515    DO  i = 0, pdims(1)-1
516       nxlf(i)   = i * nnx
517       nxrf(i)   = ( i + 1 ) * nnx - 1
518    ENDDO
519
520!
521!-- Calculate array bounds in y-direction for every PE.
522    IF ( MOD( ny+1 , pdims(2) ) /= 0 )  THEN
523       WRITE( message_string, * ) 'y-direction: gridpoint number (',ny+1,') ', &
524                           'is not an& integral divisor of the number of',      &
525                           'processors (', pdims(2),')'
526       CALL message( 'init_pegrid', 'PA0227', 1, 2, 0, 6, 0 )
527    ELSE
528       nny  = ( ny + 1 ) / pdims(2)
529    ENDIF   
530
531!
532!-- South and north array bounds
533    DO  j = 0, pdims(2)-1
534       nysf(j)   = j * nny
535       nynf(j)   = ( j + 1 ) * nny - 1
536    ENDDO
537
538!
539!-- Local array bounds of the respective PEs
540    nxl = nxlf(pcoord(1))
541    nxr = nxrf(pcoord(1))
542    nys = nysf(pcoord(2))
543    nyn = nynf(pcoord(2))
544    nzb = 0
545    nzt = nz
546    nnz = nz
547
548!
549!-- Set switches to define if the PE is situated at the border of the virtual
550!-- processor grid
551    IF ( nxl == 0 )   left_border_pe  = .TRUE.
552    IF ( nxr == nx )  right_border_pe = .TRUE.
553    IF ( nys == 0 )   south_border_pe = .TRUE.
554    IF ( nyn == ny )  north_border_pe = .TRUE.
555
556!
557!-- Calculate array bounds and gridpoint numbers for the transposed arrays
558!-- (needed in the pressure solver)
559!-- For the transposed arrays, cyclic boundaries as well as top and bottom
560!-- boundaries are omitted, because they are obstructive to the transposition
561
562!
563!-- 1. transposition  z --> x
564!-- This transposition is not neccessary in case of a 1d-decomposition along x
565    IF ( psolver == 'poisfft'  .OR.  calculate_spectra )  THEN
566
567       IF ( pdims(2) /= 1 )  THEN
568          IF ( MOD( nz , pdims(1) ) /= 0 )  THEN
569             WRITE( message_string, * ) 'transposition z --> x:',              &
570                       '& nz=',nz,' is not an integral divisior of pdims(1)=', &
571                                                                   pdims(1)
572             CALL message( 'init_pegrid', 'PA0230', 1, 2, 0, 6, 0 )
573          ENDIF
574       ENDIF
575
576       nys_x = nys
577       nyn_x = nyn
578       nny_x = nny
579       nnz_x = nz / pdims(1)
580       nzb_x = 1 + myidx * nnz_x
581       nzt_x = ( myidx + 1 ) * nnz_x
582       sendrecvcount_zx = nnx * nny * nnz_x
583
584    ENDIF
585
586
587    IF ( psolver == 'poisfft' )  THEN 
588!
589!--    2. transposition  x --> y
590       IF ( MOD( nx+1 , pdims(2) ) /= 0 )  THEN
591          WRITE( message_string, * ) 'transposition x --> y:',                 &
592                            '& nx+1=',nx+1,' is not an integral divisor of ',  &
593                            'pdims(2)=',pdims(2)
594          CALL message( 'init_pegrid', 'PA0231', 1, 2, 0, 6, 0 )
595       ENDIF
596
597       nnz_y = nnz_x
598       nzb_y = nzb_x
599       nzt_y = nzt_x
600       nnx_y = (nx+1) / pdims(2)
601       nxl_y = myidy * nnx_y
602       nxr_y = ( myidy + 1 ) * nnx_y - 1
603       sendrecvcount_xy = nnx_y * nny_x * nnz_y
604!
605!--    3. transposition  y --> z 
606!--    (ELSE:  x --> y  in case of 1D-decomposition along x)
607       nxl_z = nxl_y
608       nxr_z = nxr_y
609       nny_z = (ny+1) / pdims(1)
610       nys_z = myidx * nny_z
611       nyn_z = ( myidx + 1 ) * nny_z - 1
612       sendrecvcount_yz = nnx_y * nny_z * nnz_y
613
614       IF ( pdims(2) /= 1 )  THEN
615!
616!--       y --> z
617!--       This transposition is not neccessary in case of a 1d-decomposition
618!--       along x, except that the uptream-spline method is switched on
619          IF ( MOD( ny+1 , pdims(1) ) /= 0 )  THEN
620             WRITE( message_string, * ) 'transposition y --> z:',              &
621                               '& ny+1=',ny+1,' is not an integral divisor of',&
622                               ' pdims(1)=',pdims(1)
623             CALL message( 'init_pegrid', 'PA0232', 1, 2, 0, 6, 0 )
624          ENDIF
625
626       ELSE
627!
628!--       x --> y
629!--       This condition must be fulfilled for a 1D-decomposition along x
630          IF ( MOD( ny+1 , pdims(1) ) /= 0 )  THEN
631             WRITE( message_string, * ) 'transposition x --> y:',              &
632                               '& ny+1=',ny+1,' is not an integral divisor of',&
633                               ' pdims(1)=',pdims(1)
634             CALL message( 'init_pegrid', 'PA0233', 1, 2, 0, 6, 0 )
635          ENDIF
636
637       ENDIF
638
639    ENDIF
640
641!
642!-- Indices for direct transpositions z --> y (used for calculating spectra)
643    IF ( calculate_spectra )  THEN
644       IF ( MOD( nz, pdims(2) ) /= 0 )  THEN
645          WRITE( message_string, * ) 'direct transposition z --> y (needed ',  &
646                    'for spectra): nz=',nz,' is not an integral divisor of ',  &
647                    'pdims(2)=',pdims(2)
648          CALL message( 'init_pegrid', 'PA0234', 1, 2, 0, 6, 0 )
649       ELSE
650          nxl_yd = nxl
651          nxr_yd = nxr
652          nzb_yd = 1 + myidy * ( nz / pdims(2) )
653          nzt_yd = ( myidy + 1 ) * ( nz / pdims(2) )
654          sendrecvcount_zyd = nnx * nny * ( nz / pdims(2) )
655       ENDIF
656    ENDIF
657
658    IF ( psolver == 'poisfft'  .OR.  calculate_spectra )  THEN
659!
660!--    Indices for direct transpositions y --> x
661!--    (they are only possible in case of a 1d-decomposition along x)
662       IF ( pdims(2) == 1 )  THEN
663          nny_x = nny / pdims(1)
664          nys_x = myid * nny_x
665          nyn_x = ( myid + 1 ) * nny_x - 1
666          nzb_x = 1
667          nzt_x = nz
668          sendrecvcount_xy = nnx * nny_x * nz
669       ENDIF
670
671    ENDIF
672
673    IF ( psolver == 'poisfft' )  THEN
674!
675!--    Indices for direct transpositions x --> y
676!--    (they are only possible in case of a 1d-decomposition along y)
677       IF ( pdims(1) == 1 )  THEN
678          nnx_y = nnx / pdims(2)
679          nxl_y = myid * nnx_y
680          nxr_y = ( myid + 1 ) * nnx_y - 1
681          nzb_y = 1
682          nzt_y = nz
683          sendrecvcount_xy = nnx_y * nny * nz
684       ENDIF
685
686    ENDIF
687
688!
689!-- Arrays for storing the array bounds are needed any more
690    DEALLOCATE( nxlf , nxrf , nynf , nysf )
691
692
693!
694!-- Collect index bounds from other PEs (to be written to restart file later)
695    ALLOCATE( hor_index_bounds(4,0:numprocs-1) )
696
697    IF ( myid == 0 )  THEN
698
699       hor_index_bounds(1,0) = nxl
700       hor_index_bounds(2,0) = nxr
701       hor_index_bounds(3,0) = nys
702       hor_index_bounds(4,0) = nyn
703
704!
705!--    Receive data from all other PEs
706       DO  i = 1, numprocs-1
707          CALL MPI_RECV( ibuf, 4, MPI_INTEGER, i, MPI_ANY_TAG, comm2d, status, &
708                         ierr )
709          hor_index_bounds(:,i) = ibuf(1:4)
710       ENDDO
711
712    ELSE
713!
714!--    Send index bounds to PE0
715       ibuf(1) = nxl
716       ibuf(2) = nxr
717       ibuf(3) = nys
718       ibuf(4) = nyn
719       CALL MPI_SEND( ibuf, 4, MPI_INTEGER, 0, myid, comm2d, ierr )
720
721    ENDIF
722
723
724#if defined( __print )
725!
726!-- Control output
727    IF ( myid == 0 )  THEN
728       PRINT*, '*** processor topology ***'
729       PRINT*, ' '
730       PRINT*, 'myid   pcoord    left right  south north  idx idy   nxl: nxr',&
731               &'   nys: nyn'
732       PRINT*, '------------------------------------------------------------',&
733               &'-----------'
734       WRITE (*,1000)  0, pcoord(1), pcoord(2), pleft, pright, psouth, pnorth, &
735                       myidx, myidy, nxl, nxr, nys, nyn
7361000   FORMAT (I4,2X,'(',I3,',',I3,')',3X,I4,2X,I4,3X,I4,2X,I4,2X,I3,1X,I3, &
737               2(2X,I4,':',I4))
738
739!
740!--    Receive data from the other PEs
741       DO  i = 1,numprocs-1
742          CALL MPI_RECV( ibuf, 12, MPI_INTEGER, i, MPI_ANY_TAG, comm2d, status, &
743                         ierr )
744          WRITE (*,1000)  i, ( ibuf(j) , j = 1,12 )
745       ENDDO
746    ELSE
747
748!
749!--    Send data to PE0
750       ibuf(1) = pcoord(1); ibuf(2) = pcoord(2); ibuf(3) = pleft
751       ibuf(4) = pright; ibuf(5) = psouth; ibuf(6) = pnorth; ibuf(7) = myidx
752       ibuf(8) = myidy; ibuf(9) = nxl; ibuf(10) = nxr; ibuf(11) = nys
753       ibuf(12) = nyn
754       CALL MPI_SEND( ibuf, 12, MPI_INTEGER, 0, myid, comm2d, ierr )       
755    ENDIF
756#endif
757
758!
759!-- Determine the number of ghost point layers
760    IF ( scalar_advec   == 'ws-scheme'  .OR.                                   &
761         momentum_advec == 'ws-scheme'  .OR.  nested_run )  THEN
762       nbgp = 3
763    ELSE
764       nbgp = 1
765    ENDIF 
766
767!
768!-- Create a new MPI derived datatype for the exchange of surface (xy) data,
769!-- which is needed for coupled atmosphere-ocean runs.
770!-- First, calculate number of grid points of an xy-plane.
771    ngp_xy  = ( nxr - nxl + 1 + 2 * nbgp ) * ( nyn - nys + 1 + 2 * nbgp )
772    CALL MPI_TYPE_VECTOR( ngp_xy, 1, nzt-nzb+2, MPI_REAL, type_xy, ierr )
773    CALL MPI_TYPE_COMMIT( type_xy, ierr )
774
775    IF ( TRIM( coupling_mode ) /= 'uncoupled' .AND. .NOT. vnested )  THEN
776   
777!
778!--    Pass the number of grid points of the atmosphere model to
779!--    the ocean model and vice versa
780       IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
781
782          nx_a = nx
783          ny_a = ny
784
785          IF ( myid == 0 )  THEN
786
787             CALL MPI_SEND( nx_a, 1, MPI_INTEGER, numprocs, 1, comm_inter,  &
788                            ierr )
789             CALL MPI_SEND( ny_a, 1, MPI_INTEGER, numprocs, 2, comm_inter,  &
790                            ierr )
791             CALL MPI_SEND( pdims, 2, MPI_INTEGER, numprocs, 3, comm_inter, &
792                            ierr )
793             CALL MPI_RECV( nx_o, 1, MPI_INTEGER, numprocs, 4, comm_inter,  &
794                            status, ierr )
795             CALL MPI_RECV( ny_o, 1, MPI_INTEGER, numprocs, 5, comm_inter,  &
796                            status, ierr )
797             CALL MPI_RECV( pdims_remote, 2, MPI_INTEGER, numprocs, 6,      &
798                            comm_inter, status, ierr )
799          ENDIF
800
801          CALL MPI_BCAST( nx_o, 1, MPI_INTEGER, 0, comm2d, ierr )
802          CALL MPI_BCAST( ny_o, 1, MPI_INTEGER, 0, comm2d, ierr ) 
803          CALL MPI_BCAST( pdims_remote, 2, MPI_INTEGER, 0, comm2d, ierr )
804       
805       ELSEIF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
806
807          nx_o = nx
808          ny_o = ny 
809
810          IF ( myid == 0 ) THEN
811
812             CALL MPI_RECV( nx_a, 1, MPI_INTEGER, 0, 1, comm_inter, status, &
813                            ierr )
814             CALL MPI_RECV( ny_a, 1, MPI_INTEGER, 0, 2, comm_inter, status, &
815                            ierr )
816             CALL MPI_RECV( pdims_remote, 2, MPI_INTEGER, 0, 3, comm_inter, &
817                            status, ierr )
818             CALL MPI_SEND( nx_o, 1, MPI_INTEGER, 0, 4, comm_inter, ierr )
819             CALL MPI_SEND( ny_o, 1, MPI_INTEGER, 0, 5, comm_inter, ierr )
820             CALL MPI_SEND( pdims, 2, MPI_INTEGER, 0, 6, comm_inter, ierr )
821          ENDIF
822
823          CALL MPI_BCAST( nx_a, 1, MPI_INTEGER, 0, comm2d, ierr)
824          CALL MPI_BCAST( ny_a, 1, MPI_INTEGER, 0, comm2d, ierr) 
825          CALL MPI_BCAST( pdims_remote, 2, MPI_INTEGER, 0, comm2d, ierr) 
826
827       ENDIF
828 
829       ngp_a = ( nx_a+1 + 2 * nbgp ) * ( ny_a+1 + 2 * nbgp )
830       ngp_o = ( nx_o+1 + 2 * nbgp ) * ( ny_o+1 + 2 * nbgp )
831
832!
833!--    Determine if the horizontal grid and the number of PEs in ocean and
834!--    atmosphere is same or not
835       IF ( nx_o == nx_a  .AND.  ny_o == ny_a  .AND.  &
836            pdims(1) == pdims_remote(1) .AND. pdims(2) == pdims_remote(2) ) &
837       THEN
838          coupling_topology = 0
839       ELSE
840          coupling_topology = 1
841       ENDIF 
842
843!
844!--    Determine the target PEs for the exchange between ocean and
845!--    atmosphere (comm2d)
846       IF ( coupling_topology == 0 )  THEN
847!
848!--       In case of identical topologies, every atmosphere PE has exactly one
849!--       ocean PE counterpart and vice versa
850          IF ( TRIM( coupling_mode ) == 'atmosphere_to_ocean' ) THEN
851             target_id = myid + numprocs
852          ELSE
853             target_id = myid 
854          ENDIF
855
856       ELSE
857!
858!--       In case of nonequivalent topology in ocean and atmosphere only for
859!--       PE0 in ocean and PE0 in atmosphere a target_id is needed, since
860!--       data echxchange between ocean and atmosphere will be done only
861!--       between these PEs.   
862          IF ( myid == 0 )  THEN
863
864             IF ( TRIM( coupling_mode ) == 'atmosphere_to_ocean' )  THEN
865                target_id = numprocs 
866             ELSE
867                target_id = 0
868             ENDIF
869
870          ENDIF
871
872       ENDIF
873
874    ENDIF
875
876!
877!-- Store partner grid point co-ordinates as lists.
878!-- Create custom MPI vector datatypes for contiguous data transfer
879    IF ( vnested )  CALL vnest_init_pegrid_domain
880
881#else
882
883!
884!-- Array bounds when running on a single PE (respectively a non-parallel
885!-- machine)
886    nxl = 0
887    nxr = nx
888    nnx = nxr - nxl + 1
889    nys = 0
890    nyn = ny
891    nny = nyn - nys + 1
892    nzb = 0
893    nzt = nz
894    nnz = nz
895
896    ALLOCATE( hor_index_bounds(4,0:0) )
897    hor_index_bounds(1,0) = nxl
898    hor_index_bounds(2,0) = nxr
899    hor_index_bounds(3,0) = nys
900    hor_index_bounds(4,0) = nyn
901
902!
903!-- Array bounds for the pressure solver (in the parallel code, these bounds
904!-- are the ones for the transposed arrays)
905    nys_x = nys
906    nyn_x = nyn
907    nzb_x = nzb + 1
908    nzt_x = nzt
909
910    nxl_y = nxl
911    nxr_y = nxr
912    nzb_y = nzb + 1
913    nzt_y = nzt
914
915    nxl_z = nxl
916    nxr_z = nxr
917    nys_z = nys
918    nyn_z = nyn
919
920#endif
921
922!
923!-- Calculate number of grid levels necessary for the multigrid poisson solver
924!-- as well as the gridpoint indices on each level
925    IF ( psolver(1:9) == 'multigrid' )  THEN
926
927!
928!--    First calculate number of possible grid levels for the subdomains
929       mg_levels_x = 1
930       mg_levels_y = 1
931       mg_levels_z = 1
932
933       i = nnx
934       DO WHILE ( MOD( i, 2 ) == 0  .AND.  i /= 2 )
935          i = i / 2
936          mg_levels_x = mg_levels_x + 1
937       ENDDO
938
939       j = nny
940       DO WHILE ( MOD( j, 2 ) == 0  .AND.  j /= 2 )
941          j = j / 2
942          mg_levels_y = mg_levels_y + 1
943       ENDDO
944
945       k = nz    ! do not use nnz because it might be > nz due to transposition
946                 ! requirements
947       DO WHILE ( MOD( k, 2 ) == 0  .AND.  k /= 2 )
948          k = k / 2
949          mg_levels_z = mg_levels_z + 1
950       ENDDO
951!
952!--    The optimized MG-solver does not allow odd values for nz at the coarsest
953!--    grid level
954       IF ( TRIM( psolver ) /= 'multigrid_noopt' )  THEN
955          IF ( MOD( k, 2 ) /= 0 )  mg_levels_z = mg_levels_z - 1
956!
957!--       An odd value of nz does not work. The finest level must have an even
958!--       value.
959          IF (  mg_levels_z == 0 )  THEN
960             message_string = 'optimized multigrid method requires nz to be even'
961             CALL message( 'init_pegrid', 'PA0495', 1, 2, 0, 6, 0 )
962          ENDIF
963       ENDIF
964
965       maximum_grid_level = MIN( mg_levels_x, mg_levels_y, mg_levels_z )
966!       
967!--    Check if subdomain sizes prevents any coarsening.
968!--    This case, the maximum number of grid levels is 1, i.e. effectively
969!--    a Gauss-Seidel scheme is applied rather than a multigrid approach.
970!--    Give a warning in this case.
971       IF ( maximum_grid_level == 1  .AND.  mg_switch_to_pe0_level == -1 )  THEN
972          message_string = 'No grid coarsening possible, multigrid ' //        &
973                           'approach effectively reduces to a Gauss-Seidel ' //&
974                           'scheme.'
975         
976          CALL message( 'poismg', 'PA0648', 0, 1, 0, 6, 0 )
977       ENDIF
978
979!
980!--    Find out, if the total domain allows more levels. These additional
981!--    levels are identically processed on all PEs.
982       IF ( numprocs > 1  .AND.  mg_switch_to_pe0_level /= -1 )  THEN
983
984          IF ( mg_levels_z > MIN( mg_levels_x, mg_levels_y ) )  THEN
985
986             mg_switch_to_pe0_level_l = maximum_grid_level
987
988             mg_levels_x = 1
989             mg_levels_y = 1
990
991             i = nx+1
992             DO WHILE ( MOD( i, 2 ) == 0  .AND.  i /= 2 )
993                i = i / 2
994                mg_levels_x = mg_levels_x + 1
995             ENDDO
996
997             j = ny+1
998             DO WHILE ( MOD( j, 2 ) == 0  .AND.  j /= 2 )
999                j = j / 2
1000                mg_levels_y = mg_levels_y + 1
1001             ENDDO
1002
1003             maximum_grid_level_l = MIN( mg_levels_x, mg_levels_y, mg_levels_z )
1004
1005             IF ( maximum_grid_level_l > mg_switch_to_pe0_level_l )  THEN
1006                mg_switch_to_pe0_level_l = maximum_grid_level_l - &
1007                                           mg_switch_to_pe0_level_l + 1
1008             ELSE
1009                mg_switch_to_pe0_level_l = 0
1010             ENDIF
1011
1012          ELSE
1013
1014             mg_switch_to_pe0_level_l = 0
1015             maximum_grid_level_l = maximum_grid_level
1016
1017          ENDIF
1018
1019!
1020!--       Use switch level calculated above only if it is not pre-defined
1021!--       by user
1022          IF ( mg_switch_to_pe0_level == 0 )  THEN
1023             IF ( mg_switch_to_pe0_level_l /= 0 )  THEN
1024                mg_switch_to_pe0_level = mg_switch_to_pe0_level_l
1025                maximum_grid_level     = maximum_grid_level_l
1026             ENDIF
1027
1028          ELSE
1029!
1030!--          Check pre-defined value and reset to default, if neccessary
1031             IF ( mg_switch_to_pe0_level < mg_switch_to_pe0_level_l  .OR.      &
1032                  mg_switch_to_pe0_level >= maximum_grid_level_l )  THEN
1033                message_string = 'mg_switch_to_pe0_level ' //                  &
1034                                 'out of range and reset to 0'
1035                CALL message( 'init_pegrid', 'PA0235', 0, 1, 0, 6, 0 )
1036                mg_switch_to_pe0_level = 0
1037             ELSE
1038!
1039!--             Use the largest number of possible levels anyway and recalculate
1040!--             the switch level to this largest number of possible values
1041                maximum_grid_level = maximum_grid_level_l
1042
1043             ENDIF
1044
1045          ENDIF
1046
1047       ENDIF
1048
1049       ALLOCATE( grid_level_count(maximum_grid_level),                       &
1050                 nxl_mg(0:maximum_grid_level), nxr_mg(0:maximum_grid_level), &
1051                 nyn_mg(0:maximum_grid_level), nys_mg(0:maximum_grid_level), &
1052                 nzt_mg(0:maximum_grid_level) )
1053
1054       grid_level_count = 0
1055!
1056!--    Index zero required as dummy due to definition of arrays f2 and p2 in
1057!--    recursive subroutine next_mg_level
1058       nxl_mg(0) = 0; nxr_mg(0) = 0; nyn_mg(0) = 0; nys_mg(0) = 0; nzt_mg(0) = 0
1059
1060       nxl_l = nxl; nxr_l = nxr; nys_l = nys; nyn_l = nyn; nzt_l = nzt
1061
1062       DO  i = maximum_grid_level, 1 , -1
1063
1064          IF ( i == mg_switch_to_pe0_level )  THEN
1065#if defined( __parallel )
1066!
1067!--          Save the grid size of the subdomain at the switch level, because
1068!--          it is needed in poismg.
1069             ind(1) = nxl_l; ind(2) = nxr_l
1070             ind(3) = nys_l; ind(4) = nyn_l
1071             ind(5) = nzt_l
1072             ALLOCATE( ind_all(5*numprocs), mg_loc_ind(5,0:numprocs-1) )
1073             CALL MPI_ALLGATHER( ind, 5, MPI_INTEGER, ind_all, 5, &
1074                                 MPI_INTEGER, comm2d, ierr )
1075             DO  j = 0, numprocs-1
1076                DO  k = 1, 5
1077                   mg_loc_ind(k,j) = ind_all(k+j*5)
1078                ENDDO
1079             ENDDO
1080             DEALLOCATE( ind_all )
1081!
1082!--          Calculate the grid size of the total domain
1083             nxr_l = ( nxr_l-nxl_l+1 ) * pdims(1) - 1
1084             nxl_l = 0
1085             nyn_l = ( nyn_l-nys_l+1 ) * pdims(2) - 1
1086             nys_l = 0
1087!
1088!--          The size of this gathered array must not be larger than the
1089!--          array tend, which is used in the multigrid scheme as a temporary
1090!--          array. Therefore the subdomain size of an PE is calculated and
1091!--          the size of the gathered grid. These values are used in 
1092!--          routines pres and poismg
1093             subdomain_size = ( nxr - nxl + 2 * nbgp + 1 ) * &
1094                              ( nyn - nys + 2 * nbgp + 1 ) * ( nzt - nzb + 2 )
1095             gathered_size  = ( nxr_l - nxl_l + 3 ) * ( nyn_l - nys_l + 3 ) *  &
1096                              ( nzt_l - nzb + 2 )
1097
1098#else
1099             message_string = 'multigrid gather/scatter impossible ' //        &
1100                          'in non parallel mode'
1101             CALL message( 'init_pegrid', 'PA0237', 1, 2, 0, 6, 0 )
1102#endif
1103          ENDIF
1104
1105          nxl_mg(i) = nxl_l
1106          nxr_mg(i) = nxr_l
1107          nys_mg(i) = nys_l
1108          nyn_mg(i) = nyn_l
1109          nzt_mg(i) = nzt_l
1110
1111          nxl_l = nxl_l / 2 
1112          nxr_l = nxr_l / 2
1113          nys_l = nys_l / 2 
1114          nyn_l = nyn_l / 2 
1115          nzt_l = nzt_l / 2 
1116         
1117       ENDDO
1118
1119!
1120!--    Temporary problem: Currently calculation of maxerror in routine poismg crashes
1121!--    if grid data are collected on PE0 already on the finest grid level.
1122!--    To be solved later.
1123       IF ( maximum_grid_level == mg_switch_to_pe0_level )  THEN
1124          message_string = 'grid coarsening on subdomain level cannot be performed'
1125          CALL message( 'poismg', 'PA0236', 1, 2, 0, 6, 0 )
1126       ENDIF
1127
1128    ELSE
1129
1130       maximum_grid_level = 0
1131
1132    ENDIF
1133
1134!
1135!-- Default level 0 tells exchange_horiz that all ghost planes have to be
1136!-- exchanged. grid_level is adjusted in poismg, where only one ghost plane
1137!-- is required.
1138    grid_level = 0
1139
1140#if defined( __parallel )
1141!
1142!-- Gridpoint number for the exchange of ghost points (y-line for 2D-arrays)
1143    ngp_y  = nyn - nys + 1 + 2 * nbgp
1144
1145!
1146!-- Define new MPI derived datatypes for the exchange of ghost points in
1147!-- x- and y-direction for 2D-arrays (line)
1148    CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp, ngp_y, MPI_REAL, type_x,     &
1149                          ierr )
1150    CALL MPI_TYPE_COMMIT( type_x, ierr )
1151
1152    CALL MPI_TYPE_VECTOR( nbgp, ngp_y, ngp_y, MPI_REAL, type_y, ierr )
1153    CALL MPI_TYPE_COMMIT( type_y, ierr )
1154!
1155!-- Define new MPI derived datatypes for the exchange of ghost points in
1156!-- x- and y-direction for 2D-INTEGER arrays (line) - on normal grid.
1157!-- Define types for 32-bit and 8-bit Integer. The 8-bit Integer are only
1158!-- required on normal grid, while 32-bit Integer may be also required on
1159!-- coarser grid level in case of multigrid solver.
1160!
1161!-- 8-bit Integer
1162    CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp, ngp_y, MPI_BYTE,             &
1163                          type_x_byte, ierr )
1164    CALL MPI_TYPE_COMMIT( type_x_byte, ierr )
1165
1166    CALL MPI_TYPE_VECTOR( nbgp, ngp_y, ngp_y, MPI_BYTE,                        &
1167                          type_y_byte, ierr )
1168    CALL MPI_TYPE_COMMIT( type_y_byte, ierr )
1169!
1170!-- 32-bit Integer
1171    ALLOCATE( type_x_int(0:maximum_grid_level),                                &
1172              type_y_int(0:maximum_grid_level) )
1173             
1174    CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp, ngp_y, MPI_INTEGER,          &
1175                          type_x_int(0), ierr )
1176    CALL MPI_TYPE_COMMIT( type_x_int(0), ierr )
1177
1178    CALL MPI_TYPE_VECTOR( nbgp, ngp_y, ngp_y, MPI_INTEGER, type_y_int(0), ierr )
1179    CALL MPI_TYPE_COMMIT( type_y_int(0), ierr )
1180!
1181!-- Calculate gridpoint numbers for the exchange of ghost points along x
1182!-- (yz-plane for 3D-arrays) and define MPI derived data type(s) for the
1183!-- exchange of ghost points in y-direction (xz-plane).
1184!-- Do these calculations for the model grid and (if necessary) also
1185!-- for the coarser grid levels used in the multigrid method
1186    ALLOCATE ( ngp_xz(0:maximum_grid_level),                                   &
1187               ngp_xz_int(0:maximum_grid_level),                               &
1188               ngp_yz(0:maximum_grid_level),                                   &
1189               ngp_yz_int(0:maximum_grid_level),                               &
1190               type_xz(0:maximum_grid_level),                                  &
1191               type_xz_int(0:maximum_grid_level),                              &
1192               type_yz(0:maximum_grid_level),                                  &
1193               type_yz_int(0:maximum_grid_level) )
1194
1195    nxl_l = nxl; nxr_l = nxr; nys_l = nys; nyn_l = nyn; nzb_l = nzb; nzt_l = nzt
1196
1197!
1198!-- Discern between the model grid, which needs nbgp ghost points and
1199!-- grid levels for the multigrid scheme. In the latter case only one
1200!-- ghost point is necessary.
1201!-- First definition of MPI-datatypes for exchange of ghost layers on normal
1202!-- grid. The following loop is needed for data exchange in poismg.f90.
1203!
1204!-- Determine number of grid points of yz-layer for exchange
1205    ngp_yz(0) = (nzt - nzb + 2) * (nyn - nys + 1 + 2 * nbgp)
1206
1207!
1208!-- Define an MPI-datatype for the exchange of left/right boundaries.
1209!-- Although data are contiguous in physical memory (which does not
1210!-- necessarily require an MPI-derived datatype), the data exchange between
1211!-- left and right PE's using the MPI-derived type is 10% faster than without.
1212    CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp*(nzt-nzb+2), ngp_yz(0), &
1213                          MPI_REAL, type_xz(0), ierr )
1214    CALL MPI_TYPE_COMMIT( type_xz(0), ierr )
1215
1216    CALL MPI_TYPE_VECTOR( nbgp, ngp_yz(0), ngp_yz(0), MPI_REAL, type_yz(0), &
1217                          ierr ) 
1218    CALL MPI_TYPE_COMMIT( type_yz(0), ierr )
1219
1220!
1221!-- Define data types for exchange of 3D Integer arrays.
1222    ngp_yz_int(0) = (nzt - nzb + 2) * (nyn - nys + 1 + 2 * nbgp)
1223
1224    CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp*(nzt-nzb+2), ngp_yz_int(0),   &
1225                          MPI_INTEGER, type_xz_int(0), ierr )
1226    CALL MPI_TYPE_COMMIT( type_xz_int(0), ierr )
1227
1228    CALL MPI_TYPE_VECTOR( nbgp, ngp_yz_int(0), ngp_yz_int(0), MPI_INTEGER,     &
1229                          type_yz_int(0), ierr )
1230    CALL MPI_TYPE_COMMIT( type_yz_int(0), ierr )
1231
1232!
1233!-- Definition of MPI-datatypes for multigrid method (coarser level grids)
1234    IF ( psolver(1:9) == 'multigrid' )  THEN
1235!   
1236!--    Definition of MPI-datatyoe as above, but only 1 ghost level is used
1237       DO  i = maximum_grid_level, 1 , -1
1238!
1239!--       For 3D-exchange on different multigrid level, one ghost point for
1240!--       REAL arrays, two ghost points for INTEGER arrays
1241          ngp_xz(i) = (nzt_l - nzb_l + 2) * (nxr_l - nxl_l + 3)
1242          ngp_yz(i) = (nzt_l - nzb_l + 2) * (nyn_l - nys_l + 3)
1243
1244          ngp_xz_int(i) = (nzt_l - nzb_l + 2) * (nxr_l - nxl_l + 3)
1245          ngp_yz_int(i) = (nzt_l - nzb_l + 2) * (nyn_l - nys_l + 3)
1246!
1247!--       MPI data type for REAL arrays, for xz-layers
1248          CALL MPI_TYPE_VECTOR( nxr_l-nxl_l+3, nzt_l-nzb_l+2, ngp_yz(i),       &
1249                                MPI_REAL, type_xz(i), ierr )
1250          CALL MPI_TYPE_COMMIT( type_xz(i), ierr )
1251
1252!
1253!--       MPI data type for INTEGER arrays, for xz-layers
1254          CALL MPI_TYPE_VECTOR( nxr_l-nxl_l+3, nzt_l-nzb_l+2, ngp_yz_int(i),   &
1255                                MPI_INTEGER, type_xz_int(i), ierr )
1256          CALL MPI_TYPE_COMMIT( type_xz_int(i), ierr )
1257
1258!
1259!--       MPI data type for REAL arrays, for yz-layers
1260          CALL MPI_TYPE_VECTOR( 1, ngp_yz(i), ngp_yz(i), MPI_REAL, type_yz(i), &
1261                                ierr )
1262          CALL MPI_TYPE_COMMIT( type_yz(i), ierr )
1263!
1264!--       MPI data type for INTEGER arrays, for yz-layers
1265          CALL MPI_TYPE_VECTOR( 1, ngp_yz_int(i), ngp_yz_int(i), MPI_INTEGER,  &
1266                                type_yz_int(i), ierr )
1267          CALL MPI_TYPE_COMMIT( type_yz_int(i), ierr )
1268
1269
1270!--       For 2D-exchange of INTEGER arrays on coarser grid level, where 2 ghost
1271!--       points need to be exchanged. Only required for 32-bit Integer arrays.
1272          CALL MPI_TYPE_VECTOR( nxr_l-nxl_l+5, 2, nyn_l-nys_l+5, MPI_INTEGER,  &
1273                                type_x_int(i), ierr )
1274          CALL MPI_TYPE_COMMIT( type_x_int(i), ierr )
1275
1276
1277          CALL MPI_TYPE_VECTOR( 2, nyn_l-nys_l+5, nyn_l-nys_l+5, MPI_INTEGER,  &
1278                                type_y_int(i), ierr )
1279          CALL MPI_TYPE_COMMIT( type_y_int(i), ierr )
1280
1281          nxl_l = nxl_l / 2
1282          nxr_l = nxr_l / 2
1283          nys_l = nys_l / 2
1284          nyn_l = nyn_l / 2
1285          nzt_l = nzt_l / 2
1286
1287       ENDDO
1288
1289    ENDIF
1290
1291#endif
1292
1293#if defined( __parallel )
1294!
1295!-- Setting of flags for inflow/outflow/nesting conditions.
1296    IF ( pleft == MPI_PROC_NULL )  THEN
1297       IF ( bc_lr == 'dirichlet/radiation'  .OR.  bc_lr == 'nested'  .OR.      &
1298            bc_lr == 'nesting_offline' )  THEN
1299          bc_dirichlet_l  = .TRUE.
1300       ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
1301          bc_radiation_l = .TRUE.
1302       ENDIF
1303    ENDIF
1304 
1305    IF ( pright == MPI_PROC_NULL )  THEN
1306       IF ( bc_lr == 'dirichlet/radiation' )  THEN
1307          bc_radiation_r = .TRUE.
1308       ELSEIF ( bc_lr == 'radiation/dirichlet'  .OR.  bc_lr == 'nested'  .OR.  &
1309                bc_lr == 'nesting_offline' )  THEN
1310          bc_dirichlet_r  = .TRUE.
1311       ENDIF
1312    ENDIF
1313
1314    IF ( psouth == MPI_PROC_NULL )  THEN
1315       IF ( bc_ns == 'dirichlet/radiation' )  THEN
1316          bc_radiation_s = .TRUE.
1317       ELSEIF ( bc_ns == 'radiation/dirichlet'  .OR.  bc_ns == 'nested'  .OR.  &
1318                bc_ns == 'nesting_offline' )  THEN
1319          bc_dirichlet_s  = .TRUE.
1320       ENDIF
1321    ENDIF
1322
1323    IF ( pnorth == MPI_PROC_NULL )  THEN
1324       IF ( bc_ns == 'dirichlet/radiation'  .OR.  bc_ns == 'nested'  .OR.      &
1325            bc_ns == 'nesting_offline' )  THEN
1326          bc_dirichlet_n  = .TRUE.
1327       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
1328          bc_radiation_n = .TRUE.
1329       ENDIF
1330    ENDIF
1331!
1332!-- In case of synthetic turbulence geneartor determine ids.
1333!-- Please note, if no forcing or nesting is applied, the generator is applied
1334!-- only at the left lateral boundary.
1335    IF ( use_syn_turb_gen )  THEN
1336       IF ( bc_dirichlet_l )  THEN
1337          id_stg_left_l = myidx
1338       ELSE
1339          id_stg_left_l = 0
1340       ENDIF
1341       IF ( bc_dirichlet_r )  THEN
1342          id_stg_right_l = myidx
1343       ELSE
1344          id_stg_right_l = 0
1345       ENDIF
1346       IF ( bc_dirichlet_s )  THEN
1347          id_stg_south_l = myidy
1348       ELSE
1349          id_stg_south_l = 0
1350       ENDIF
1351       IF ( bc_dirichlet_n )  THEN
1352          id_stg_north_l = myidy
1353       ELSE
1354          id_stg_north_l = 0
1355       ENDIF
1356
1357       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1358       CALL MPI_ALLREDUCE( id_stg_left_l, id_stg_left,   1, MPI_INTEGER,       &
1359                           MPI_SUM, comm1dx, ierr )
1360
1361       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1362       CALL MPI_ALLREDUCE( id_stg_right_l, id_stg_right, 1, MPI_INTEGER,       &
1363                           MPI_SUM, comm1dx, ierr )
1364
1365       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1366       CALL MPI_ALLREDUCE( id_stg_south_l, id_stg_south, 1, MPI_INTEGER,       &
1367                           MPI_SUM, comm1dy, ierr )
1368
1369       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1370       CALL MPI_ALLREDUCE( id_stg_north_l, id_stg_north, 1, MPI_INTEGER,       &
1371                           MPI_SUM, comm1dy, ierr )
1372
1373    ENDIF 
1374 
1375!
1376!-- Broadcast the id of the inflow PE
1377    IF ( bc_dirichlet_l )  THEN
1378       id_inflow_l = myidx
1379    ELSE
1380       id_inflow_l = 0
1381    ENDIF
1382    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1383    CALL MPI_ALLREDUCE( id_inflow_l, id_inflow, 1, MPI_INTEGER, MPI_SUM, &
1384                        comm1dx, ierr )
1385
1386!
1387!-- Broadcast the id of the recycling plane
1388!-- WARNING: needs to be adjusted in case of inflows other than from left side!
1389    IF ( turbulent_inflow ) THEN
1390   
1391       IF ( NINT( recycling_width / dx, KIND=idp ) >= nxl  .AND.                                   &
1392            NINT( recycling_width / dx, KIND=idp ) <= nxr )  THEN
1393          id_recycling_l = myidx
1394       ELSE
1395          id_recycling_l = 0
1396       ENDIF
1397       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1398       CALL MPI_ALLREDUCE( id_recycling_l, id_recycling, 1, MPI_INTEGER, MPI_SUM, &
1399                           comm1dx, ierr )
1400                           
1401    ENDIF
1402
1403!
1404!-- Broadcast the id of the outflow PE and outflow-source plane
1405    IF ( turbulent_outflow )  THEN
1406
1407       IF ( bc_radiation_r )  THEN
1408          id_outflow_l = myidx
1409       ELSE
1410          id_outflow_l = 0
1411       ENDIF
1412       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1413       CALL MPI_ALLREDUCE( id_outflow_l, id_outflow, 1, MPI_INTEGER, MPI_SUM, &
1414                           comm1dx, ierr )
1415
1416       IF ( NINT( outflow_source_plane / dx ) >= nxl  .AND. &
1417            NINT( outflow_source_plane / dx ) <= nxr )  THEN
1418          id_outflow_source_l = myidx
1419       ELSE
1420          id_outflow_source_l = 0
1421       ENDIF
1422       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1423       CALL MPI_ALLREDUCE( id_outflow_source_l, id_outflow_source, 1, &
1424                           MPI_INTEGER, MPI_SUM, comm1dx, ierr )
1425
1426    ENDIF
1427
1428    CALL location_message( 'creating virtual PE grids + MPI derived data types', 'finished' )
1429
1430#else
1431    IF ( bc_lr == 'dirichlet/radiation' )  THEN
1432       bc_dirichlet_l = .TRUE.
1433       bc_radiation_r = .TRUE.
1434    ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
1435       bc_radiation_l = .TRUE.
1436       bc_dirichlet_r = .TRUE.
1437    ENDIF
1438
1439    IF ( bc_ns == 'dirichlet/radiation' )  THEN
1440       bc_dirichlet_n = .TRUE.
1441       bc_radiation_s = .TRUE.
1442    ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
1443       bc_radiation_n = .TRUE.
1444       bc_dirichlet_s = .TRUE.
1445    ENDIF
1446#endif
1447
1448!
1449!-- At the inflow or outflow, u or v, respectively, have to be calculated for
1450!-- one more grid point.
1451    IF ( bc_dirichlet_l  .OR.  bc_radiation_l )  THEN
1452       nxlu = nxl + 1
1453    ELSE
1454       nxlu = nxl
1455    ENDIF
1456    IF ( bc_dirichlet_s  .OR.  bc_radiation_s )  THEN
1457       nysv = nys + 1
1458    ELSE
1459       nysv = nys
1460    ENDIF
1461
1462!
1463!-- Allocate wall flag arrays used in the multigrid solver
1464    IF ( psolver(1:9) == 'multigrid' )  THEN
1465
1466       DO  i = maximum_grid_level, 1, -1
1467
1468           SELECT CASE ( i )
1469
1470              CASE ( 1 )
1471                 ALLOCATE( wall_flags_1(nzb:nzt_mg(i)+1,         &
1472                                        nys_mg(i)-1:nyn_mg(i)+1, &
1473                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1474
1475              CASE ( 2 )
1476                 ALLOCATE( wall_flags_2(nzb:nzt_mg(i)+1,         &
1477                                        nys_mg(i)-1:nyn_mg(i)+1, &
1478                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1479
1480              CASE ( 3 )
1481                 ALLOCATE( wall_flags_3(nzb:nzt_mg(i)+1,         &
1482                                        nys_mg(i)-1:nyn_mg(i)+1, &
1483                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1484
1485              CASE ( 4 )
1486                 ALLOCATE( wall_flags_4(nzb:nzt_mg(i)+1,         &
1487                                        nys_mg(i)-1:nyn_mg(i)+1, &
1488                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1489
1490              CASE ( 5 )
1491                 ALLOCATE( wall_flags_5(nzb:nzt_mg(i)+1,         &
1492                                        nys_mg(i)-1:nyn_mg(i)+1, &
1493                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1494
1495              CASE ( 6 )
1496                 ALLOCATE( wall_flags_6(nzb:nzt_mg(i)+1,         &
1497                                        nys_mg(i)-1:nyn_mg(i)+1, &
1498                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1499
1500              CASE ( 7 )
1501                 ALLOCATE( wall_flags_7(nzb:nzt_mg(i)+1,         &
1502                                        nys_mg(i)-1:nyn_mg(i)+1, &
1503                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1504
1505              CASE ( 8 )
1506                 ALLOCATE( wall_flags_8(nzb:nzt_mg(i)+1,         &
1507                                        nys_mg(i)-1:nyn_mg(i)+1, &
1508                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1509
1510              CASE ( 9 )
1511                 ALLOCATE( wall_flags_9(nzb:nzt_mg(i)+1,         &
1512                                        nys_mg(i)-1:nyn_mg(i)+1, &
1513                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1514
1515              CASE ( 10 )
1516                 ALLOCATE( wall_flags_10(nzb:nzt_mg(i)+1,        &
1517                                        nys_mg(i)-1:nyn_mg(i)+1, &
1518                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1519
1520              CASE DEFAULT
1521                 message_string = 'more than 10 multigrid levels'
1522                 CALL message( 'init_pegrid', 'PA0238', 1, 2, 0, 6, 0 )
1523
1524          END SELECT
1525
1526       ENDDO
1527
1528    ENDIF
1529
1530 END SUBROUTINE init_pegrid
Note: See TracBrowser for help on using the repository browser.