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

Last change on this file since 1734 was 1683, checked in by knoop, 9 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 42.3 KB
Line 
1!> @file init_pegrid.f90
2!--------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms
6! of the GNU General Public License as published by the Free Software Foundation,
7! either version 3 of the License, or (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
10! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
11! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with
14! PALM. If not, see <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2014 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------!
18!
19! Current revisions:
20! ------------------
21!
22!
23! Former revisions:
24! -----------------
25! $Id: init_pegrid.f90 1683 2015-10-07 23:57:51Z raasch $
26!
27! 1682 2015-10-07 23:56:08Z knoop
28! Code annotations made doxygen readable
29!
30! 1677 2015-10-02 13:25:23Z boeske
31! New MPI-data types for exchange of 3D integer arrays.
32!
33! 1575 2015-03-27 09:56:27Z raasch
34! adjustments for psolver-queries, calculation of ngp_xz added
35!
36! 1565 2015-03-09 20:59:31Z suehring
37! Refine if-clause for setting nbgp.
38!
39! 1557 2015-03-05 16:43:04Z suehring
40! Adjustment for monotonic limiter
41!
42! 1468 2014-09-24 14:06:57Z maronga
43! Adapted for use on up to 6-digit processor cores
44!
45! 1435 2014-07-21 10:37:02Z keck
46! bugfix: added missing parameter coupling_mode_remote to ONLY-attribute
47!
48! 1402 2014-05-09 14:25:13Z raasch
49! location messages modified
50!
51! 1384 2014-05-02 14:31:06Z raasch
52! location messages added
53!
54! 1353 2014-04-08 15:21:23Z heinze
55! REAL constants provided with KIND-attribute
56!
57! 1322 2014-03-20 16:38:49Z raasch
58! REAL functions provided with KIND-attribute
59!
60! 1320 2014-03-20 08:40:49Z raasch
61! ONLY-attribute added to USE-statements,
62! kind-parameters added to all INTEGER and REAL declaration statements,
63! kinds are defined in new module kinds,
64! revision history before 2012 removed,
65! comment fields (!:) to be used for variable explanations added to
66! all variable declaration statements
67!
68! 1304 2014-03-12 10:29:42Z raasch
69! bugfix: single core MPI runs missed some settings of transpose indices
70!
71! 1212 2013-08-15 08:46:27Z raasch
72! error message for poisfft_hybrid removed
73!
74! 1159 2013-05-21 11:58:22Z fricke
75! dirichlet/neumann and neumann/dirichlet removed
76!
77! 1139 2013-04-18 07:25:03Z raasch
78! bugfix for calculating the id of the PE carrying the recycling plane
79!
80! 1111 2013-03-08 23:54:10Z raasch
81! initialization of poisfft moved to module poisfft
82!
83! 1092 2013-02-02 11:24:22Z raasch
84! unused variables removed
85!
86! 1056 2012-11-16 15:28:04Z raasch
87! Indices for arrays n.._mg start from zero due to definition of arrays f2 and
88! p2 as automatic arrays in recursive subroutine next_mg_level
89!
90! 1041 2012-11-06 02:36:29Z raasch
91! a 2d virtual processor topology is used by default for all machines
92!
93! 1036 2012-10-22 13:43:42Z raasch
94! code put under GPL (PALM 3.9)
95!
96! 1003 2012-09-14 14:35:53Z raasch
97! subdomains must have identical size (grid matching = "match" removed)
98!
99! 1001 2012-09-13 14:08:46Z raasch
100! all actions concerning upstream-spline-method removed
101!
102! 978 2012-08-09 08:28:32Z fricke
103! dirichlet/neumann and neumann/dirichlet added
104! nxlu and nysv are also calculated for inflow boundary
105!
106! 809 2012-01-30 13:32:58Z maronga
107! Bugfix: replaced .AND. and .NOT. with && and ! in the preprocessor directives
108!
109! 807 2012-01-25 11:53:51Z maronga
110! New cpp directive "__check" implemented which is used by check_namelist_files
111!
112! Revision 1.1  1997/07/24 11:15:09  raasch
113! Initial revision
114!
115!
116! Description:
117! ------------
118!> Determination of the virtual processor topology (if not prescribed by the
119!> user)and computation of the grid point number and array bounds of the local
120!> domains.
121!------------------------------------------------------------------------------!
122 SUBROUTINE init_pegrid
123 
124
125    USE control_parameters,                                                    &
126        ONLY:  bc_lr, bc_ns, coupling_mode, coupling_mode_remote,              &
127               coupling_topology, dt_dosp, gathered_size, grid_level,          &
128               grid_level_count, host, inflow_l, inflow_n, inflow_r, inflow_s, &
129               io_blocks, io_group, maximum_grid_level,                        &
130               maximum_parallel_io_streams, message_string,                    &
131               mg_switch_to_pe0_level, momentum_advec, neutral, psolver,       &
132               outflow_l, outflow_n, outflow_r, outflow_s, recycling_width,    &
133               scalar_advec, subdomain_size 
134
135    USE grid_variables,                                                        &
136        ONLY:  dx
137       
138    USE indices,                                                               &
139        ONLY:  mg_loc_ind, nbgp, nnx, nny, nnz, nx, nx_a, nx_o, nxl, nxl_mg,   &
140               nxlu, nxr, nxr_mg, ny, ny_a, ny_o, nyn, nyn_mg, nys, nys_mg,    &
141               nysv, nz, nzb, nzt, nzt_mg, wall_flags_1, wall_flags_2,         &
142               wall_flags_3, wall_flags_4, wall_flags_5, wall_flags_6,         &
143               wall_flags_7, wall_flags_8, wall_flags_9, wall_flags_10
144
145    USE kinds
146     
147    USE pegrid
148 
149    USE transpose_indices,                                                     &
150        ONLY:  nxl_y, nxl_yd, nxl_z, nxr_y, nxr_yd, nxr_z, nyn_x, nyn_z, nys_x,&
151               nys_z, nzb_x, nzb_y, nzb_yd, nzt_x, nzt_yd, nzt_y
152
153    IMPLICIT NONE
154
155    INTEGER(iwp) ::  i                        !<
156    INTEGER(iwp) ::  id_inflow_l              !<
157    INTEGER(iwp) ::  id_recycling_l           !<
158    INTEGER(iwp) ::  ind(5)                   !<
159    INTEGER(iwp) ::  j                        !<
160    INTEGER(iwp) ::  k                        !<
161    INTEGER(iwp) ::  maximum_grid_level_l     !<
162    INTEGER(iwp) ::  mg_levels_x              !<
163    INTEGER(iwp) ::  mg_levels_y              !<
164    INTEGER(iwp) ::  mg_levels_z              !<
165    INTEGER(iwp) ::  mg_switch_to_pe0_level_l !<
166    INTEGER(iwp) ::  nnx_y                    !<
167    INTEGER(iwp) ::  nnx_z                    !<
168    INTEGER(iwp) ::  nny_x                    !<
169    INTEGER(iwp) ::  nny_z                    !<
170    INTEGER(iwp) ::  nnz_x                    !<
171    INTEGER(iwp) ::  nnz_y                    !<
172    INTEGER(iwp) ::  numproc_sqr              !<
173    INTEGER(iwp) ::  nxl_l                    !<
174    INTEGER(iwp) ::  nxr_l                    !<
175    INTEGER(iwp) ::  nyn_l                    !<
176    INTEGER(iwp) ::  nys_l                    !<
177    INTEGER(iwp) ::  nzb_l                    !<
178    INTEGER(iwp) ::  nzt_l                    !<
179    INTEGER(iwp) ::  omp_get_num_threads      !<
180
181    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  ind_all !<
182    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nxlf    !<
183    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nxrf    !<
184    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nynf    !<
185    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nysf    !<
186
187    INTEGER(iwp), DIMENSION(2) :: pdims_remote          !<
188
189#if defined( __mpi2 )
190    LOGICAL ::  found                                   !<
191#endif
192
193!
194!-- Get the number of OpenMP threads
195    !$OMP PARALLEL
196#if defined( __intel_openmp_bug )
197    threads_per_task = omp_get_num_threads()
198#else
199!$  threads_per_task = omp_get_num_threads()
200#endif
201    !$OMP END PARALLEL
202
203
204#if defined( __parallel )
205
206    CALL location_message( 'creating virtual PE grids + MPI derived data types', &
207                           .FALSE. )
208!
209!-- Determine the processor topology or check it, if prescribed by the user
210    IF ( npex == -1  .AND.  npey == -1 )  THEN
211
212!
213!--    Automatic determination of the topology
214       numproc_sqr = SQRT( REAL( numprocs, KIND=wp ) )
215       pdims(1)    = MAX( numproc_sqr , 1 )
216       DO  WHILE ( MOD( numprocs , pdims(1) ) /= 0 )
217          pdims(1) = pdims(1) - 1
218       ENDDO
219       pdims(2) = numprocs / pdims(1)
220
221    ELSEIF ( npex /= -1  .AND.  npey /= -1 )  THEN
222
223!
224!--    Prescribed by user. Number of processors on the prescribed topology
225!--    must be equal to the number of PEs available to the job
226       IF ( ( npex * npey ) /= numprocs )  THEN
227          WRITE( message_string, * ) 'number of PEs of the prescribed ',      & 
228                 'topology (', npex*npey,') does not match & the number of ', & 
229                 'PEs available to the job (', numprocs, ')'
230          CALL message( 'init_pegrid', 'PA0221', 1, 2, 0, 6, 0 )
231       ENDIF
232       pdims(1) = npex
233       pdims(2) = npey
234
235    ELSE
236!
237!--    If the processor topology is prescribed by the user, the number of
238!--    PEs must be given in both directions
239       message_string = 'if the processor topology is prescribed by the, ' //  &
240                   ' user& both values of "npex" and "npey" must be given ' // &
241                   'in the &NAMELIST-parameter file'
242       CALL message( 'init_pegrid', 'PA0222', 1, 2, 0, 6, 0 )
243
244    ENDIF
245
246!
247!-- For communication speedup, set barriers in front of collective
248!-- communications by default on SGI-type systems
249    IF ( host(3:5) == 'sgi' )  collective_wait = .TRUE.
250
251!
252!-- If necessary, set horizontal boundary conditions to non-cyclic
253    IF ( bc_lr /= 'cyclic' )  cyclic(1) = .FALSE.
254    IF ( bc_ns /= 'cyclic' )  cyclic(2) = .FALSE.
255
256
257#if ! defined( __check)
258!
259!-- Create the virtual processor grid
260    CALL MPI_CART_CREATE( comm_palm, ndim, pdims, cyclic, reorder, &
261                          comm2d, ierr )
262    CALL MPI_COMM_RANK( comm2d, myid, ierr )
263    WRITE (myid_char,'(''_'',I6.6)')  myid
264
265    CALL MPI_CART_COORDS( comm2d, myid, ndim, pcoord, ierr )
266    CALL MPI_CART_SHIFT( comm2d, 0, 1, pleft, pright, ierr )
267    CALL MPI_CART_SHIFT( comm2d, 1, 1, psouth, pnorth, ierr )
268
269!
270!-- Determine sub-topologies for transpositions
271!-- Transposition from z to x:
272    remain_dims(1) = .TRUE.
273    remain_dims(2) = .FALSE.
274    CALL MPI_CART_SUB( comm2d, remain_dims, comm1dx, ierr )
275    CALL MPI_COMM_RANK( comm1dx, myidx, ierr )
276!
277!-- Transposition from x to y
278    remain_dims(1) = .FALSE.
279    remain_dims(2) = .TRUE.
280    CALL MPI_CART_SUB( comm2d, remain_dims, comm1dy, ierr )
281    CALL MPI_COMM_RANK( comm1dy, myidy, ierr )
282
283#endif
284
285!
286!-- Calculate array bounds along x-direction for every PE.
287    ALLOCATE( nxlf(0:pdims(1)-1), nxrf(0:pdims(1)-1), nynf(0:pdims(2)-1), &
288              nysf(0:pdims(2)-1) )
289
290    IF ( MOD( nx+1 , pdims(1) ) /= 0 )  THEN
291       WRITE( message_string, * ) 'x-direction: gridpoint number (',nx+1,') ',&
292                               'is not an& integral divisor of the number ',  &
293                               'processors (', pdims(1),')'
294       CALL message( 'init_pegrid', 'PA0225', 1, 2, 0, 6, 0 )
295    ELSE
296       nnx  = ( nx + 1 ) / pdims(1)
297       IF ( nnx*pdims(1) - ( nx + 1) > nnx )  THEN
298          WRITE( message_string, * ) 'x-direction: nx does not match the',    & 
299                       'requirements given by the number of PEs &used',       &
300                       '& please use nx = ', nx - ( pdims(1) - ( nnx*pdims(1) &
301                                   - ( nx + 1 ) ) ), ' instead of nx =', nx
302          CALL message( 'init_pegrid', 'PA0226', 1, 2, 0, 6, 0 )
303       ENDIF
304    ENDIF   
305
306!
307!-- Left and right array bounds, number of gridpoints
308    DO  i = 0, pdims(1)-1
309       nxlf(i)   = i * nnx
310       nxrf(i)   = ( i + 1 ) * nnx - 1
311    ENDDO
312
313!
314!-- Calculate array bounds in y-direction for every PE.
315    IF ( MOD( ny+1 , pdims(2) ) /= 0 )  THEN
316       WRITE( message_string, * ) 'y-direction: gridpoint number (',ny+1,') ', &
317                           'is not an& integral divisor of the number of',     &
318                           'processors (', pdims(2),')'
319       CALL message( 'init_pegrid', 'PA0227', 1, 2, 0, 6, 0 )
320    ELSE
321       nny  = ( ny + 1 ) / pdims(2)
322       IF ( nny*pdims(2) - ( ny + 1) > nny )  THEN
323          WRITE( message_string, * ) 'y-direction: ny does not match the',    &
324                       'requirements given by the number of PEs &used ',      &
325                       '& please use ny = ', ny - ( pdims(2) - ( nnx*pdims(2) &
326                                     - ( ny + 1 ) ) ), ' instead of ny =', ny
327          CALL message( 'init_pegrid', 'PA0228', 1, 2, 0, 6, 0 )
328       ENDIF
329    ENDIF   
330
331!
332!-- South and north array bounds
333    DO  j = 0, pdims(2)-1
334       nysf(j)   = j * nny
335       nynf(j)   = ( j + 1 ) * nny - 1
336    ENDDO
337
338!
339!-- Local array bounds of the respective PEs
340    nxl = nxlf(pcoord(1))
341    nxr = nxrf(pcoord(1))
342    nys = nysf(pcoord(2))
343    nyn = nynf(pcoord(2))
344    nzb = 0
345    nzt = nz
346    nnz = nz
347
348!
349!-- Set switches to define if the PE is situated at the border of the virtual
350!-- processor grid
351    IF ( nxl == 0 )   left_border_pe  = .TRUE.
352    IF ( nxr == nx )  right_border_pe = .TRUE.
353    IF ( nys == 0 )   south_border_pe = .TRUE.
354    IF ( nyn == ny )  north_border_pe = .TRUE.
355
356!
357!-- Calculate array bounds and gridpoint numbers for the transposed arrays
358!-- (needed in the pressure solver)
359!-- For the transposed arrays, cyclic boundaries as well as top and bottom
360!-- boundaries are omitted, because they are obstructive to the transposition
361
362!
363!-- 1. transposition  z --> x
364!-- This transposition is not neccessary in case of a 1d-decomposition along x
365    nys_x = nys
366    nyn_x = nyn
367    nny_x = nny
368    nnz_x = nz / pdims(1)
369    nzb_x = 1 + myidx * nnz_x
370    nzt_x = ( myidx + 1 ) * nnz_x
371    sendrecvcount_zx = nnx * nny * nnz_x
372
373    IF ( pdims(2) /= 1 )  THEN
374       IF ( MOD( nz , pdims(1) ) /= 0 )  THEN
375          WRITE( message_string, * ) 'transposition z --> x:',                &
376                       '&nz=',nz,' is not an integral divisior of pdims(1)=', &
377                                                                   pdims(1)
378          CALL message( 'init_pegrid', 'PA0230', 1, 2, 0, 6, 0 )
379       ENDIF
380    ENDIF
381
382!
383!-- 2. transposition  x --> y
384    nnz_y = nnz_x
385    nzb_y = nzb_x
386    nzt_y = nzt_x
387    IF ( MOD( nx+1 , pdims(2) ) /= 0 )  THEN
388       WRITE( message_string, * ) 'transposition x --> y:',                &
389                         '&nx+1=',nx+1,' is not an integral divisor of ',&
390                         'pdims(2)=',pdims(2)
391       CALL message( 'init_pegrid', 'PA0231', 1, 2, 0, 6, 0 )
392    ENDIF
393    nnx_y = (nx+1) / pdims(2)
394    nxl_y = myidy * nnx_y
395    nxr_y = ( myidy + 1 ) * nnx_y - 1
396    sendrecvcount_xy = nnx_y * nny_x * nnz_y
397
398!
399!-- 3. transposition  y --> z  (ELSE:  x --> y  in case of 1D-decomposition
400!-- along x)
401    nnx_z = nnx_y
402    nxl_z = nxl_y
403    nxr_z = nxr_y
404    nny_z = (ny+1) / pdims(1)
405    nys_z = myidx * nny_z
406    nyn_z = ( myidx + 1 ) * nny_z - 1
407    sendrecvcount_yz = nnx_y * nny_z * nnz_y
408
409    IF ( pdims(2) /= 1 )  THEN
410!
411!--    y --> z
412!--    This transposition is not neccessary in case of a 1d-decomposition
413!--    along x, except that the uptream-spline method is switched on
414       IF ( MOD( ny+1 , pdims(1) ) /= 0 )  THEN
415          WRITE( message_string, * ) 'transposition y --> z:',            &
416                            '& ny+1=',ny+1,' is not an integral divisor of ',&
417                            'pdims(1)=',pdims(1)
418          CALL message( 'init_pegrid', 'PA0232', 1, 2, 0, 6, 0 )
419       ENDIF
420
421    ELSE
422!
423!--    x --> y. This condition must be fulfilled for a 1D-decomposition along x
424       IF ( MOD( ny+1 , pdims(1) ) /= 0 )  THEN
425          WRITE( message_string, * ) 'transposition x --> y:',               &
426                            '& ny+1=',ny+1,' is not an integral divisor of ',&
427                            'pdims(1)=',pdims(1)
428          CALL message( 'init_pegrid', 'PA0233', 1, 2, 0, 6, 0 )
429       ENDIF
430
431    ENDIF
432
433!
434!-- Indices for direct transpositions z --> y (used for calculating spectra)
435    IF ( dt_dosp /= 9999999.9_wp )  THEN
436       IF ( MOD( nz, pdims(2) ) /= 0 )  THEN
437          WRITE( message_string, * ) 'direct transposition z --> y (needed ', &
438                    'for spectra):& nz=',nz,' is not an integral divisor of ',&
439                    'pdims(2)=',pdims(2)
440          CALL message( 'init_pegrid', 'PA0234', 1, 2, 0, 6, 0 )
441       ELSE
442          nxl_yd = nxl
443          nxr_yd = nxr
444          nzb_yd = 1 + myidy * ( nz / pdims(2) )
445          nzt_yd = ( myidy + 1 ) * ( nz / pdims(2) )
446          sendrecvcount_zyd = nnx * nny * ( nz / pdims(2) )
447       ENDIF
448    ENDIF
449
450!
451!-- Indices for direct transpositions y --> x (they are only possible in case
452!-- of a 1d-decomposition along x)
453    IF ( pdims(2) == 1 )  THEN
454       nny_x = nny / pdims(1)
455       nys_x = myid * nny_x
456       nyn_x = ( myid + 1 ) * nny_x - 1
457       nzb_x = 1
458       nzt_x = nz
459       sendrecvcount_xy = nnx * nny_x * nz
460    ENDIF
461
462!
463!-- Indices for direct transpositions x --> y (they are only possible in case
464!-- of a 1d-decomposition along y)
465    IF ( pdims(1) == 1 )  THEN
466       nnx_y = nnx / pdims(2)
467       nxl_y = myid * nnx_y
468       nxr_y = ( myid + 1 ) * nnx_y - 1
469       nzb_y = 1
470       nzt_y = nz
471       sendrecvcount_xy = nnx_y * nny * nz
472    ENDIF
473
474!
475!-- Arrays for storing the array bounds are needed any more
476    DEALLOCATE( nxlf , nxrf , nynf , nysf )
477
478
479#if ! defined( __check)
480!
481!-- Collect index bounds from other PEs (to be written to restart file later)
482    ALLOCATE( hor_index_bounds(4,0:numprocs-1) )
483
484    IF ( myid == 0 )  THEN
485
486       hor_index_bounds(1,0) = nxl
487       hor_index_bounds(2,0) = nxr
488       hor_index_bounds(3,0) = nys
489       hor_index_bounds(4,0) = nyn
490
491!
492!--    Receive data from all other PEs
493       DO  i = 1, numprocs-1
494          CALL MPI_RECV( ibuf, 4, MPI_INTEGER, i, MPI_ANY_TAG, comm2d, status, &
495                         ierr )
496          hor_index_bounds(:,i) = ibuf(1:4)
497       ENDDO
498
499    ELSE
500!
501!--    Send index bounds to PE0
502       ibuf(1) = nxl
503       ibuf(2) = nxr
504       ibuf(3) = nys
505       ibuf(4) = nyn
506       CALL MPI_SEND( ibuf, 4, MPI_INTEGER, 0, myid, comm2d, ierr )
507
508    ENDIF
509
510#endif
511
512#if defined( __print )
513!
514!-- Control output
515    IF ( myid == 0 )  THEN
516       PRINT*, '*** processor topology ***'
517       PRINT*, ' '
518       PRINT*, 'myid   pcoord    left right  south north  idx idy   nxl: nxr',&
519               &'   nys: nyn'
520       PRINT*, '------------------------------------------------------------',&
521               &'-----------'
522       WRITE (*,1000)  0, pcoord(1), pcoord(2), pleft, pright, psouth, pnorth, &
523                       myidx, myidy, nxl, nxr, nys, nyn
5241000   FORMAT (I4,2X,'(',I3,',',I3,')',3X,I4,2X,I4,3X,I4,2X,I4,2X,I3,1X,I3, &
525               2(2X,I4,':',I4))
526
527!
528!--    Receive data from the other PEs
529       DO  i = 1,numprocs-1
530          CALL MPI_RECV( ibuf, 12, MPI_INTEGER, i, MPI_ANY_TAG, comm2d, status, &
531                         ierr )
532          WRITE (*,1000)  i, ( ibuf(j) , j = 1,12 )
533       ENDDO
534    ELSE
535
536!
537!--    Send data to PE0
538       ibuf(1) = pcoord(1); ibuf(2) = pcoord(2); ibuf(3) = pleft
539       ibuf(4) = pright; ibuf(5) = psouth; ibuf(6) = pnorth; ibuf(7) = myidx
540       ibuf(8) = myidy; ibuf(9) = nxl; ibuf(10) = nxr; ibuf(11) = nys
541       ibuf(12) = nyn
542       CALL MPI_SEND( ibuf, 12, MPI_INTEGER, 0, myid, comm2d, ierr )       
543    ENDIF
544#endif
545
546#if defined( __parallel ) && ! defined( __check)
547#if defined( __mpi2 )
548!
549!-- In case of coupled runs, get the port name on PE0 of the atmosphere model
550!-- and pass it to PE0 of the ocean model
551    IF ( myid == 0 )  THEN
552
553       IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
554
555          CALL MPI_OPEN_PORT( MPI_INFO_NULL, port_name, ierr )
556
557          CALL MPI_PUBLISH_NAME( 'palm_coupler', MPI_INFO_NULL, port_name, &
558                                 ierr )
559
560!
561!--       Write a flag file for the ocean model and the other atmosphere
562!--       processes.
563!--       There seems to be a bug in MPICH2 which causes hanging processes
564!--       in case that execution of LOOKUP_NAME is continued too early
565!--       (i.e. before the port has been created)
566          OPEN( 90, FILE='COUPLING_PORT_OPENED', FORM='FORMATTED' )
567          WRITE ( 90, '(''TRUE'')' )
568          CLOSE ( 90 )
569
570       ELSEIF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
571
572!
573!--       Continue only if the atmosphere model has created the port.
574!--       There seems to be a bug in MPICH2 which causes hanging processes
575!--       in case that execution of LOOKUP_NAME is continued too early
576!--       (i.e. before the port has been created)
577          INQUIRE( FILE='COUPLING_PORT_OPENED', EXIST=found )
578          DO WHILE ( .NOT. found )
579             INQUIRE( FILE='COUPLING_PORT_OPENED', EXIST=found )
580          ENDDO
581
582          CALL MPI_LOOKUP_NAME( 'palm_coupler', MPI_INFO_NULL, port_name, ierr )
583
584       ENDIF
585
586    ENDIF
587
588!
589!-- In case of coupled runs, establish the connection between the atmosphere
590!-- and the ocean model and define the intercommunicator (comm_inter)
591    CALL MPI_BARRIER( comm2d, ierr )
592    IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
593
594       CALL MPI_COMM_ACCEPT( port_name, MPI_INFO_NULL, 0, MPI_COMM_WORLD, &
595                             comm_inter, ierr )
596       coupling_mode_remote = 'ocean_to_atmosphere'
597
598    ELSEIF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
599
600       CALL MPI_COMM_CONNECT( port_name, MPI_INFO_NULL, 0, MPI_COMM_WORLD, &
601                              comm_inter, ierr )
602       coupling_mode_remote = 'atmosphere_to_ocean'
603
604    ENDIF
605#endif
606
607!
608!-- Determine the number of ghost point layers
609    IF ( ( scalar_advec == 'ws-scheme' .AND. .NOT. neutral ) .OR.             &
610         scalar_advec == 'ws-scheme-mono' .OR.                                &
611         momentum_advec == 'ws-scheme' )  THEN
612       nbgp = 3
613    ELSE
614       nbgp = 1
615    ENDIF 
616
617!
618!-- Create a new MPI derived datatype for the exchange of surface (xy) data,
619!-- which is needed for coupled atmosphere-ocean runs.
620!-- First, calculate number of grid points of an xy-plane.
621    ngp_xy  = ( nxr - nxl + 1 + 2 * nbgp ) * ( nyn - nys + 1 + 2 * nbgp )
622    CALL MPI_TYPE_VECTOR( ngp_xy, 1, nzt-nzb+2, MPI_REAL, type_xy, ierr )
623    CALL MPI_TYPE_COMMIT( type_xy, ierr )
624
625    IF ( TRIM( coupling_mode ) /= 'uncoupled' )  THEN
626   
627!
628!--    Pass the number of grid points of the atmosphere model to
629!--    the ocean model and vice versa
630       IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
631
632          nx_a = nx
633          ny_a = ny
634
635          IF ( myid == 0 )  THEN
636
637             CALL MPI_SEND( nx_a, 1, MPI_INTEGER, numprocs, 1, comm_inter,  &
638                            ierr )
639             CALL MPI_SEND( ny_a, 1, MPI_INTEGER, numprocs, 2, comm_inter,  &
640                            ierr )
641             CALL MPI_SEND( pdims, 2, MPI_INTEGER, numprocs, 3, comm_inter, &
642                            ierr )
643             CALL MPI_RECV( nx_o, 1, MPI_INTEGER, numprocs, 4, comm_inter,  &
644                            status, ierr )
645             CALL MPI_RECV( ny_o, 1, MPI_INTEGER, numprocs, 5, comm_inter,  &
646                            status, ierr )
647             CALL MPI_RECV( pdims_remote, 2, MPI_INTEGER, numprocs, 6,      &
648                            comm_inter, status, ierr )
649          ENDIF
650
651          CALL MPI_BCAST( nx_o, 1, MPI_INTEGER, 0, comm2d, ierr )
652          CALL MPI_BCAST( ny_o, 1, MPI_INTEGER, 0, comm2d, ierr ) 
653          CALL MPI_BCAST( pdims_remote, 2, MPI_INTEGER, 0, comm2d, ierr )
654       
655       ELSEIF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
656
657          nx_o = nx
658          ny_o = ny 
659
660          IF ( myid == 0 ) THEN
661
662             CALL MPI_RECV( nx_a, 1, MPI_INTEGER, 0, 1, comm_inter, status, &
663                            ierr )
664             CALL MPI_RECV( ny_a, 1, MPI_INTEGER, 0, 2, comm_inter, status, &
665                            ierr )
666             CALL MPI_RECV( pdims_remote, 2, MPI_INTEGER, 0, 3, comm_inter, &
667                            status, ierr )
668             CALL MPI_SEND( nx_o, 1, MPI_INTEGER, 0, 4, comm_inter, ierr )
669             CALL MPI_SEND( ny_o, 1, MPI_INTEGER, 0, 5, comm_inter, ierr )
670             CALL MPI_SEND( pdims, 2, MPI_INTEGER, 0, 6, comm_inter, ierr )
671          ENDIF
672
673          CALL MPI_BCAST( nx_a, 1, MPI_INTEGER, 0, comm2d, ierr)
674          CALL MPI_BCAST( ny_a, 1, MPI_INTEGER, 0, comm2d, ierr) 
675          CALL MPI_BCAST( pdims_remote, 2, MPI_INTEGER, 0, comm2d, ierr) 
676
677       ENDIF
678 
679       ngp_a = ( nx_a+1 + 2 * nbgp ) * ( ny_a+1 + 2 * nbgp )
680       ngp_o = ( nx_o+1 + 2 * nbgp ) * ( ny_o+1 + 2 * nbgp )
681
682!
683!--    Determine if the horizontal grid and the number of PEs in ocean and
684!--    atmosphere is same or not
685       IF ( nx_o == nx_a  .AND.  ny_o == ny_a  .AND.  &
686            pdims(1) == pdims_remote(1) .AND. pdims(2) == pdims_remote(2) ) &
687       THEN
688          coupling_topology = 0
689       ELSE
690          coupling_topology = 1
691       ENDIF 
692
693!
694!--    Determine the target PEs for the exchange between ocean and
695!--    atmosphere (comm2d)
696       IF ( coupling_topology == 0 )  THEN
697!
698!--       In case of identical topologies, every atmosphere PE has exactly one
699!--       ocean PE counterpart and vice versa
700          IF ( TRIM( coupling_mode ) == 'atmosphere_to_ocean' ) THEN
701             target_id = myid + numprocs
702          ELSE
703             target_id = myid 
704          ENDIF
705
706       ELSE
707!
708!--       In case of nonequivalent topology in ocean and atmosphere only for
709!--       PE0 in ocean and PE0 in atmosphere a target_id is needed, since
710!--       data echxchange between ocean and atmosphere will be done only
711!--       between these PEs.   
712          IF ( myid == 0 )  THEN
713
714             IF ( TRIM( coupling_mode ) == 'atmosphere_to_ocean' )  THEN
715                target_id = numprocs 
716             ELSE
717                target_id = 0
718             ENDIF
719
720          ENDIF
721
722       ENDIF
723
724    ENDIF
725
726
727#endif
728
729#else
730
731!
732!-- Array bounds when running on a single PE (respectively a non-parallel
733!-- machine)
734    nxl = 0
735    nxr = nx
736    nnx = nxr - nxl + 1
737    nys = 0
738    nyn = ny
739    nny = nyn - nys + 1
740    nzb = 0
741    nzt = nz
742    nnz = nz
743
744    ALLOCATE( hor_index_bounds(4,0:0) )
745    hor_index_bounds(1,0) = nxl
746    hor_index_bounds(2,0) = nxr
747    hor_index_bounds(3,0) = nys
748    hor_index_bounds(4,0) = nyn
749
750!
751!-- Array bounds for the pressure solver (in the parallel code, these bounds
752!-- are the ones for the transposed arrays)
753    nys_x = nys
754    nyn_x = nyn
755    nzb_x = nzb + 1
756    nzt_x = nzt
757
758    nxl_y = nxl
759    nxr_y = nxr
760    nzb_y = nzb + 1
761    nzt_y = nzt
762
763    nxl_z = nxl
764    nxr_z = nxr
765    nys_z = nys
766    nyn_z = nyn
767
768#endif
769
770!
771!-- Calculate number of grid levels necessary for the multigrid poisson solver
772!-- as well as the gridpoint indices on each level
773    IF ( psolver(1:9) == 'multigrid' )  THEN
774
775!
776!--    First calculate number of possible grid levels for the subdomains
777       mg_levels_x = 1
778       mg_levels_y = 1
779       mg_levels_z = 1
780
781       i = nnx
782       DO WHILE ( MOD( i, 2 ) == 0  .AND.  i /= 2 )
783          i = i / 2
784          mg_levels_x = mg_levels_x + 1
785       ENDDO
786
787       j = nny
788       DO WHILE ( MOD( j, 2 ) == 0  .AND.  j /= 2 )
789          j = j / 2
790          mg_levels_y = mg_levels_y + 1
791       ENDDO
792
793       k = nz    ! do not use nnz because it might be > nz due to transposition
794                 ! requirements
795       DO WHILE ( MOD( k, 2 ) == 0  .AND.  k /= 2 )
796          k = k / 2
797          mg_levels_z = mg_levels_z + 1
798       ENDDO
799
800       maximum_grid_level = MIN( mg_levels_x, mg_levels_y, mg_levels_z )
801
802!
803!--    Find out, if the total domain allows more levels. These additional
804!--    levels are identically processed on all PEs.
805       IF ( numprocs > 1  .AND.  mg_switch_to_pe0_level /= -1 )  THEN
806
807          IF ( mg_levels_z > MIN( mg_levels_x, mg_levels_y ) )  THEN
808
809             mg_switch_to_pe0_level_l = maximum_grid_level
810
811             mg_levels_x = 1
812             mg_levels_y = 1
813
814             i = nx+1
815             DO WHILE ( MOD( i, 2 ) == 0  .AND.  i /= 2 )
816                i = i / 2
817                mg_levels_x = mg_levels_x + 1
818             ENDDO
819
820             j = ny+1
821             DO WHILE ( MOD( j, 2 ) == 0  .AND.  j /= 2 )
822                j = j / 2
823                mg_levels_y = mg_levels_y + 1
824             ENDDO
825
826             maximum_grid_level_l = MIN( mg_levels_x, mg_levels_y, mg_levels_z )
827
828             IF ( maximum_grid_level_l > mg_switch_to_pe0_level_l )  THEN
829                mg_switch_to_pe0_level_l = maximum_grid_level_l - &
830                                           mg_switch_to_pe0_level_l + 1
831             ELSE
832                mg_switch_to_pe0_level_l = 0
833             ENDIF
834
835          ELSE
836             mg_switch_to_pe0_level_l = 0
837             maximum_grid_level_l = maximum_grid_level
838
839          ENDIF
840
841!
842!--       Use switch level calculated above only if it is not pre-defined
843!--       by user
844          IF ( mg_switch_to_pe0_level == 0 )  THEN
845             IF ( mg_switch_to_pe0_level_l /= 0 )  THEN
846                mg_switch_to_pe0_level = mg_switch_to_pe0_level_l
847                maximum_grid_level     = maximum_grid_level_l
848             ENDIF
849
850          ELSE
851!
852!--          Check pre-defined value and reset to default, if neccessary
853             IF ( mg_switch_to_pe0_level < mg_switch_to_pe0_level_l  .OR.  &
854                  mg_switch_to_pe0_level >= maximum_grid_level_l )  THEN
855                message_string = 'mg_switch_to_pe0_level ' // &
856                                 'out of range and reset to default (=0)'
857                CALL message( 'init_pegrid', 'PA0235', 0, 1, 0, 6, 0 )
858                mg_switch_to_pe0_level = 0
859             ELSE
860!
861!--             Use the largest number of possible levels anyway and recalculate
862!--             the switch level to this largest number of possible values
863                maximum_grid_level = maximum_grid_level_l
864
865             ENDIF
866
867          ENDIF
868
869       ENDIF
870
871       ALLOCATE( grid_level_count(maximum_grid_level),                       &
872                 nxl_mg(0:maximum_grid_level), nxr_mg(0:maximum_grid_level), &
873                 nyn_mg(0:maximum_grid_level), nys_mg(0:maximum_grid_level), &
874                 nzt_mg(0:maximum_grid_level) )
875
876       grid_level_count = 0
877!
878!--    Index zero required as dummy due to definition of arrays f2 and p2 in
879!--    recursive subroutine next_mg_level
880       nxl_mg(0) = 0; nxr_mg(0) = 0; nyn_mg(0) = 0; nys_mg(0) = 0; nzt_mg(0) = 0
881
882       nxl_l = nxl; nxr_l = nxr; nys_l = nys; nyn_l = nyn; nzt_l = nzt
883
884       DO  i = maximum_grid_level, 1 , -1
885
886          IF ( i == mg_switch_to_pe0_level )  THEN
887#if defined( __parallel ) && ! defined( __check )
888!
889!--          Save the grid size of the subdomain at the switch level, because
890!--          it is needed in poismg.
891             ind(1) = nxl_l; ind(2) = nxr_l
892             ind(3) = nys_l; ind(4) = nyn_l
893             ind(5) = nzt_l
894             ALLOCATE( ind_all(5*numprocs), mg_loc_ind(5,0:numprocs-1) )
895             CALL MPI_ALLGATHER( ind, 5, MPI_INTEGER, ind_all, 5, &
896                                 MPI_INTEGER, comm2d, ierr )
897             DO  j = 0, numprocs-1
898                DO  k = 1, 5
899                   mg_loc_ind(k,j) = ind_all(k+j*5)
900                ENDDO
901             ENDDO
902             DEALLOCATE( ind_all )
903!
904!--          Calculate the grid size of the total domain
905             nxr_l = ( nxr_l-nxl_l+1 ) * pdims(1) - 1
906             nxl_l = 0
907             nyn_l = ( nyn_l-nys_l+1 ) * pdims(2) - 1
908             nys_l = 0
909!
910!--          The size of this gathered array must not be larger than the
911!--          array tend, which is used in the multigrid scheme as a temporary
912!--          array. Therefore the subdomain size of an PE is calculated and
913!--          the size of the gathered grid. These values are used in 
914!--          routines pres and poismg
915             subdomain_size = ( nxr - nxl + 2 * nbgp + 1 ) * &
916                              ( nyn - nys + 2 * nbgp + 1 ) * ( nzt - nzb + 2 )
917             gathered_size  = ( nxr_l - nxl_l + 3 ) * ( nyn_l - nys_l + 3 ) * &
918                              ( nzt_l - nzb + 2 )
919
920#elif ! defined ( __parallel )
921             message_string = 'multigrid gather/scatter impossible ' // &
922                          'in non parallel mode'
923             CALL message( 'init_pegrid', 'PA0237', 1, 2, 0, 6, 0 )
924#endif
925          ENDIF
926
927          nxl_mg(i) = nxl_l
928          nxr_mg(i) = nxr_l
929          nys_mg(i) = nys_l
930          nyn_mg(i) = nyn_l
931          nzt_mg(i) = nzt_l
932
933          nxl_l = nxl_l / 2 
934          nxr_l = nxr_l / 2
935          nys_l = nys_l / 2 
936          nyn_l = nyn_l / 2 
937          nzt_l = nzt_l / 2 
938
939       ENDDO
940
941!
942!--    Temporary problem: Currently calculation of maxerror iin routine poismg crashes
943!--    if grid data are collected on PE0 already on the finest grid level.
944!--    To be solved later.
945       IF ( maximum_grid_level == mg_switch_to_pe0_level )  THEN
946          message_string = 'grid coarsening on subdomain level cannot be performed'
947          CALL message( 'poismg', 'PA0236', 1, 2, 0, 6, 0 )
948       ENDIF
949
950    ELSE
951
952       maximum_grid_level = 0
953
954    ENDIF
955
956!
957!-- Default level 0 tells exchange_horiz that all ghost planes have to be
958!-- exchanged. grid_level is adjusted in poismg, where only one ghost plane
959!-- is required.
960    grid_level = 0
961
962#if defined( __parallel ) && ! defined ( __check )
963!
964!-- Gridpoint number for the exchange of ghost points (y-line for 2D-arrays)
965    ngp_y  = nyn - nys + 1 + 2 * nbgp
966
967!
968!-- Define new MPI derived datatypes for the exchange of ghost points in
969!-- x- and y-direction for 2D-arrays (line)
970    CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp, ngp_y, MPI_REAL, type_x, &
971                          ierr )
972    CALL MPI_TYPE_COMMIT( type_x, ierr )
973    CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp, ngp_y, MPI_INTEGER, &
974                          type_x_int, ierr )
975    CALL MPI_TYPE_COMMIT( type_x_int, ierr )
976
977    CALL MPI_TYPE_VECTOR( nbgp, ngp_y, ngp_y, MPI_REAL, type_y, ierr )
978    CALL MPI_TYPE_COMMIT( type_y, ierr )
979    CALL MPI_TYPE_VECTOR( nbgp, ngp_y, ngp_y, MPI_INTEGER, type_y_int, ierr )
980    CALL MPI_TYPE_COMMIT( type_y_int, ierr )
981
982
983!
984!-- Calculate gridpoint numbers for the exchange of ghost points along x
985!-- (yz-plane for 3D-arrays) and define MPI derived data type(s) for the
986!-- exchange of ghost points in y-direction (xz-plane).
987!-- Do these calculations for the model grid and (if necessary) also
988!-- for the coarser grid levels used in the multigrid method
989    ALLOCATE ( ngp_xz(0:maximum_grid_level), ngp_yz(0:maximum_grid_level),     &
990               type_xz(0:maximum_grid_level), type_yz(0:maximum_grid_level) )
991
992    nxl_l = nxl; nxr_l = nxr; nys_l = nys; nyn_l = nyn; nzb_l = nzb; nzt_l = nzt
993
994!
995!-- Discern between the model grid, which needs nbgp ghost points and
996!-- grid levels for the multigrid scheme. In the latter case only one
997!-- ghost point is necessary.
998!-- First definition of MPI-datatypes for exchange of ghost layers on normal
999!-- grid. The following loop is needed for data exchange in poismg.f90.
1000!
1001!-- Determine number of grid points of yz-layer for exchange
1002    ngp_yz(0) = (nzt - nzb + 2) * (nyn - nys + 1 + 2 * nbgp)
1003
1004!
1005!-- Define an MPI-datatype for the exchange of left/right boundaries.
1006!-- Although data are contiguous in physical memory (which does not
1007!-- necessarily require an MPI-derived datatype), the data exchange between
1008!-- left and right PE's using the MPI-derived type is 10% faster than without.
1009    CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp*(nzt-nzb+2), ngp_yz(0), &
1010                          MPI_REAL, type_xz(0), ierr )
1011    CALL MPI_TYPE_COMMIT( type_xz(0), ierr )
1012
1013    CALL MPI_TYPE_VECTOR( nbgp, ngp_yz(0), ngp_yz(0), MPI_REAL, type_yz(0), &
1014                          ierr ) 
1015    CALL MPI_TYPE_COMMIT( type_yz(0), ierr )
1016
1017!
1018!-- Definition of MPI-datatypes for multigrid method (coarser level grids)
1019    IF ( psolver(1:9) == 'multigrid' )  THEN
1020!   
1021!--    Definition of MPI-datatyoe as above, but only 1 ghost level is used
1022       DO  i = maximum_grid_level, 1 , -1
1023
1024          ngp_xz(i) = (nzt_l - nzb_l + 2) * (nxr_l - nxl_l + 3)
1025          ngp_yz(i) = (nzt_l - nzb_l + 2) * (nyn_l - nys_l + 3)
1026
1027          CALL MPI_TYPE_VECTOR( nxr_l-nxl_l+3, nzt_l-nzb_l+2, ngp_yz(i), &
1028                                MPI_REAL, type_xz(i), ierr )
1029          CALL MPI_TYPE_COMMIT( type_xz(i), ierr )
1030
1031          CALL MPI_TYPE_VECTOR( 1, ngp_yz(i), ngp_yz(i), MPI_REAL, type_yz(i), &
1032                                ierr )
1033          CALL MPI_TYPE_COMMIT( type_yz(i), ierr )
1034
1035          nxl_l = nxl_l / 2
1036          nxr_l = nxr_l / 2
1037          nys_l = nys_l / 2
1038          nyn_l = nyn_l / 2
1039          nzt_l = nzt_l / 2
1040
1041       ENDDO
1042
1043    ENDIF
1044!
1045!-- Define data types for exchange of 3D Integer arrays.
1046    ngp_yz_int = (nzt - nzb + 2) * (nyn - nys + 1 + 2 * nbgp)
1047
1048    CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp*(nzt-nzb+2), ngp_yz_int, &
1049                          MPI_INTEGER, type_xz_int, ierr )
1050    CALL MPI_TYPE_COMMIT( type_xz_int, ierr )
1051
1052    CALL MPI_TYPE_VECTOR( nbgp, ngp_yz_int, ngp_yz_int, MPI_INTEGER, type_yz_int, &
1053                          ierr )
1054    CALL MPI_TYPE_COMMIT( type_yz_int, ierr )
1055
1056#endif
1057
1058#if defined( __parallel ) && ! defined ( __check )
1059!
1060!-- Setting of flags for inflow/outflow conditions in case of non-cyclic
1061!-- horizontal boundary conditions.
1062    IF ( pleft == MPI_PROC_NULL )  THEN
1063       IF ( bc_lr == 'dirichlet/radiation' )  THEN
1064          inflow_l  = .TRUE.
1065       ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
1066          outflow_l = .TRUE.
1067       ENDIF
1068    ENDIF
1069
1070    IF ( pright == MPI_PROC_NULL )  THEN
1071       IF ( bc_lr == 'dirichlet/radiation' )  THEN
1072          outflow_r = .TRUE.
1073       ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
1074          inflow_r  = .TRUE.
1075       ENDIF
1076    ENDIF
1077
1078    IF ( psouth == MPI_PROC_NULL )  THEN
1079       IF ( bc_ns == 'dirichlet/radiation' )  THEN
1080          outflow_s = .TRUE.
1081       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
1082          inflow_s  = .TRUE.
1083       ENDIF
1084    ENDIF
1085
1086    IF ( pnorth == MPI_PROC_NULL )  THEN
1087       IF ( bc_ns == 'dirichlet/radiation' )  THEN
1088          inflow_n  = .TRUE.
1089       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
1090          outflow_n = .TRUE.
1091       ENDIF
1092    ENDIF
1093
1094!
1095!-- Broadcast the id of the inflow PE
1096    IF ( inflow_l )  THEN
1097       id_inflow_l = myidx
1098    ELSE
1099       id_inflow_l = 0
1100    ENDIF
1101    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1102    CALL MPI_ALLREDUCE( id_inflow_l, id_inflow, 1, MPI_INTEGER, MPI_SUM, &
1103                        comm1dx, ierr )
1104
1105!
1106!-- Broadcast the id of the recycling plane
1107!-- WARNING: needs to be adjusted in case of inflows other than from left side!
1108    IF ( NINT( recycling_width / dx ) >= nxl  .AND. &
1109         NINT( recycling_width / dx ) <= nxr )  THEN
1110       id_recycling_l = myidx
1111    ELSE
1112       id_recycling_l = 0
1113    ENDIF
1114    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1115    CALL MPI_ALLREDUCE( id_recycling_l, id_recycling, 1, MPI_INTEGER, MPI_SUM, &
1116                        comm1dx, ierr )
1117
1118    CALL location_message( 'finished', .TRUE. )
1119
1120#elif ! defined ( __parallel )
1121    IF ( bc_lr == 'dirichlet/radiation' )  THEN
1122       inflow_l  = .TRUE.
1123       outflow_r = .TRUE.
1124    ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
1125       outflow_l = .TRUE.
1126       inflow_r  = .TRUE.
1127    ENDIF
1128
1129    IF ( bc_ns == 'dirichlet/radiation' )  THEN
1130       inflow_n  = .TRUE.
1131       outflow_s = .TRUE.
1132    ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
1133       outflow_n = .TRUE.
1134       inflow_s  = .TRUE.
1135    ENDIF
1136#endif
1137
1138!
1139!-- At the inflow or outflow, u or v, respectively, have to be calculated for
1140!-- one more grid point.
1141    IF ( inflow_l .OR. outflow_l )  THEN
1142       nxlu = nxl + 1
1143    ELSE
1144       nxlu = nxl
1145    ENDIF
1146    IF ( inflow_s .OR. outflow_s )  THEN
1147       nysv = nys + 1
1148    ELSE
1149       nysv = nys
1150    ENDIF
1151
1152!
1153!-- Allocate wall flag arrays used in the multigrid solver
1154    IF ( psolver(1:9) == 'multigrid' )  THEN
1155
1156       DO  i = maximum_grid_level, 1, -1
1157
1158           SELECT CASE ( i )
1159
1160              CASE ( 1 )
1161                 ALLOCATE( wall_flags_1(nzb:nzt_mg(i)+1,         &
1162                                        nys_mg(i)-1:nyn_mg(i)+1, &
1163                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1164
1165              CASE ( 2 )
1166                 ALLOCATE( wall_flags_2(nzb:nzt_mg(i)+1,         &
1167                                        nys_mg(i)-1:nyn_mg(i)+1, &
1168                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1169
1170              CASE ( 3 )
1171                 ALLOCATE( wall_flags_3(nzb:nzt_mg(i)+1,         &
1172                                        nys_mg(i)-1:nyn_mg(i)+1, &
1173                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1174
1175              CASE ( 4 )
1176                 ALLOCATE( wall_flags_4(nzb:nzt_mg(i)+1,         &
1177                                        nys_mg(i)-1:nyn_mg(i)+1, &
1178                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1179
1180              CASE ( 5 )
1181                 ALLOCATE( wall_flags_5(nzb:nzt_mg(i)+1,         &
1182                                        nys_mg(i)-1:nyn_mg(i)+1, &
1183                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1184
1185              CASE ( 6 )
1186                 ALLOCATE( wall_flags_6(nzb:nzt_mg(i)+1,         &
1187                                        nys_mg(i)-1:nyn_mg(i)+1, &
1188                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1189
1190              CASE ( 7 )
1191                 ALLOCATE( wall_flags_7(nzb:nzt_mg(i)+1,         &
1192                                        nys_mg(i)-1:nyn_mg(i)+1, &
1193                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1194
1195              CASE ( 8 )
1196                 ALLOCATE( wall_flags_8(nzb:nzt_mg(i)+1,         &
1197                                        nys_mg(i)-1:nyn_mg(i)+1, &
1198                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1199
1200              CASE ( 9 )
1201                 ALLOCATE( wall_flags_9(nzb:nzt_mg(i)+1,         &
1202                                        nys_mg(i)-1:nyn_mg(i)+1, &
1203                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1204
1205              CASE ( 10 )
1206                 ALLOCATE( wall_flags_10(nzb:nzt_mg(i)+1,        &
1207                                        nys_mg(i)-1:nyn_mg(i)+1, &
1208                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1209
1210              CASE DEFAULT
1211                 message_string = 'more than 10 multigrid levels'
1212                 CALL message( 'init_pegrid', 'PA0238', 1, 2, 0, 6, 0 )
1213
1214          END SELECT
1215
1216       ENDDO
1217
1218    ENDIF
1219
1220!
1221!-- Calculate the number of groups into which parallel I/O is split.
1222!-- The default for files which are opened by all PEs (or where each
1223!-- PE opens his own independent file) is, that all PEs are doing input/output
1224!-- in parallel at the same time. This might cause performance or even more
1225!-- severe problems depending on the configuration of the underlying file
1226!-- system.
1227!-- First, set the default:
1228    IF ( maximum_parallel_io_streams == -1  .OR. &
1229         maximum_parallel_io_streams > numprocs )  THEN
1230       maximum_parallel_io_streams = numprocs
1231    ENDIF
1232
1233!
1234!-- Now calculate the number of io_blocks and the io_group to which the
1235!-- respective PE belongs. I/O of the groups is done in serial, but in parallel
1236!-- for all PEs belonging to the same group. A preliminary setting with myid
1237!-- based on MPI_COMM_WORLD has been done in parin.
1238    io_blocks = numprocs / maximum_parallel_io_streams
1239    io_group  = MOD( myid+1, io_blocks )
1240   
1241
1242 END SUBROUTINE init_pegrid
Note: See TracBrowser for help on using the repository browser.