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

Last change on this file since 1140 was 1140, checked in by raasch, 8 years ago

last commit documented

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