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

Last change on this file since 181 was 181, checked in by raasch, 13 years ago

bugfixes + adjustments for SGI ICE system

  • Property svn:keywords set to Id
File size: 33.2 KB
Line 
1 SUBROUTINE init_pegrid
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6! ATTENTION: nnz_x undefined problem still has to be solved!!!!!!!!
7! nz is used instead nnz for calculating mg-levels
8! Collect on PE0 horizontal index bounds from all other PEs,
9! broadcast the id of the inflow PE (using the respective communicator)
10! TEST OUTPUT (TO BE REMOVED) logging mpi2 ierr values
11!
12! Former revisions:
13! -----------------
14! $Id: init_pegrid.f90 181 2008-07-30 07:07:47Z raasch $
15!
16! 114 2007-10-10 00:03:15Z raasch
17! Allocation of wall flag arrays for multigrid solver
18!
19! 108 2007-08-24 15:10:38Z letzel
20! Intercommunicator (comm_inter) and derived data type (type_xy) for
21! coupled model runs created, assign coupling_mode_remote,
22! indices nxlu and nysv are calculated (needed for non-cyclic boundary
23! conditions)
24!
25! 82 2007-04-16 15:40:52Z raasch
26! Cpp-directive lcmuk changed to intel_openmp_bug, setting of host on lcmuk by
27! cpp-directive removed
28!
29! 75 2007-03-22 09:54:05Z raasch
30! uxrp, vynp eliminated,
31! dirichlet/neumann changed to dirichlet/radiation, etc.,
32! poisfft_init is only called if fft-solver is switched on
33!
34! RCS Log replace by Id keyword, revision history cleaned up
35!
36! Revision 1.28  2006/04/26 13:23:32  raasch
37! lcmuk does not understand the !$ comment so a cpp-directive is required
38!
39! Revision 1.1  1997/07/24 11:15:09  raasch
40! Initial revision
41!
42!
43! Description:
44! ------------
45! Determination of the virtual processor topology (if not prescribed by the
46! user)and computation of the grid point number and array bounds of the local
47! domains.
48!------------------------------------------------------------------------------!
49
50    USE control_parameters
51    USE fft_xy
52    USE grid_variables
53    USE indices
54    USE pegrid
55    USE poisfft_mod
56    USE poisfft_hybrid_mod
57    USE statistics
58    USE transpose_indices
59
60
61    IMPLICIT NONE
62
63    INTEGER ::  gathered_size, i, id_inflow_l, id_recycling_l, ind(5), j, k, &
64                maximum_grid_level_l, mg_switch_to_pe0_level_l, mg_levels_x, &
65                mg_levels_y, mg_levels_z, nnx_y, nnx_z, nny_x, nny_z, nnz_x, &
66                nnz_y, numproc_sqr, nx_total, nxl_l, nxr_l, nyn_l, nys_l,    &
67                nzb_l, nzt_l, omp_get_num_threads, subdomain_size
68
69    INTEGER, DIMENSION(:), ALLOCATABLE ::  ind_all, nxlf, nxrf, nynf, nysf
70
71    LOGICAL ::  found
72
73!
74!-- Get the number of OpenMP threads
75    !$OMP PARALLEL
76#if defined( __intel_openmp_bug )
77    threads_per_task = omp_get_num_threads()
78#else
79!$  threads_per_task = omp_get_num_threads()
80#endif
81    !$OMP END PARALLEL
82
83
84#if defined( __parallel )
85!
86!-- Determine the processor topology or check it, if prescribed by the user
87    IF ( npex == -1  .AND.  npey == -1 )  THEN
88
89!
90!--    Automatic determination of the topology
91!--    The default on SMP- and cluster-hosts is a 1d-decomposition along x
92       IF ( host(1:3) == 'ibm'  .OR.  host(1:3) == 'nec'  .OR. &
93            host(1:2) == 'lc'   .OR.  host(1:3) == 'dec' )  THEN
94
95          pdims(1) = numprocs
96          pdims(2) = 1
97
98       ELSE
99
100          numproc_sqr = SQRT( REAL( numprocs ) )
101          pdims(1)    = MAX( numproc_sqr , 1 )
102          DO  WHILE ( MOD( numprocs , pdims(1) ) /= 0 )
103             pdims(1) = pdims(1) - 1
104          ENDDO
105          pdims(2) = numprocs / pdims(1)
106
107       ENDIF
108
109    ELSEIF ( npex /= -1  .AND.  npey /= -1 )  THEN
110
111!
112!--    Prescribed by user. Number of processors on the prescribed topology
113!--    must be equal to the number of PEs available to the job
114       IF ( ( npex * npey ) /= numprocs )  THEN
115          PRINT*, '+++ init_pegrid:'
116          PRINT*, '    number of PEs of the prescribed topology (', npex*npey, &
117                      ') does not match the number of PEs available to the ',  &
118                      'job (', numprocs, ')'
119          CALL local_stop
120       ENDIF
121       pdims(1) = npex
122       pdims(2) = npey
123
124    ELSE
125!
126!--    If the processor topology is prescribed by the user, the number of
127!--    PEs must be given in both directions
128       PRINT*, '+++ init_pegrid:'
129       PRINT*, '    if the processor topology is prescribed by the user, ',   &
130                    'both values of "npex" and "npey" must be given in the ', &
131                    'NAMELIST-parameter file'
132       CALL local_stop
133
134    ENDIF
135
136!
137!-- The hybrid solver can only be used in case of a 1d-decomposition along x
138    IF ( pdims(2) /= 1  .AND.  psolver == 'poisfft_hybrid' )  THEN
139       IF ( myid == 0 )  THEN
140          PRINT*, '*** init_pegrid: psolver = "poisfft_hybrid" can only be'
141          PRINT*, '                 used in case of a 1d-decomposition along x'
142       ENDIF
143    ENDIF
144
145!
146!-- If necessary, set horizontal boundary conditions to non-cyclic
147    IF ( bc_lr /= 'cyclic' )  cyclic(1) = .FALSE.
148    IF ( bc_ns /= 'cyclic' )  cyclic(2) = .FALSE.
149
150!
151!-- Create the virtual processor grid
152    CALL MPI_CART_CREATE( comm_palm, ndim, pdims, cyclic, reorder, &
153                          comm2d, ierr )
154    CALL MPI_COMM_RANK( comm2d, myid, ierr )
155    WRITE (myid_char,'(''_'',I4.4)')  myid
156
157    CALL MPI_CART_COORDS( comm2d, myid, ndim, pcoord, ierr )
158    CALL MPI_CART_SHIFT( comm2d, 0, 1, pleft, pright, ierr )
159    CALL MPI_CART_SHIFT( comm2d, 1, 1, psouth, pnorth, ierr )
160
161!
162!-- Determine sub-topologies for transpositions
163!-- Transposition from z to x:
164    remain_dims(1) = .TRUE.
165    remain_dims(2) = .FALSE.
166    CALL MPI_CART_SUB( comm2d, remain_dims, comm1dx, ierr )
167    CALL MPI_COMM_RANK( comm1dx, myidx, ierr )
168!
169!-- Transposition from x to y
170    remain_dims(1) = .FALSE.
171    remain_dims(2) = .TRUE.
172    CALL MPI_CART_SUB( comm2d, remain_dims, comm1dy, ierr )
173    CALL MPI_COMM_RANK( comm1dy, myidy, ierr )
174
175
176!
177!-- Find a grid (used for array d) which will match the transposition demands
178    IF ( grid_matching == 'strict' )  THEN
179
180       nxa = nx;  nya = ny;  nza = nz
181
182    ELSE
183
184       found = .FALSE.
185   xn: DO  nxa = nx, 2*nx
186!
187!--       Meet conditions for nx
188          IF ( MOD( nxa+1, pdims(1) ) /= 0 .OR. &
189               MOD( nxa+1, pdims(2) ) /= 0 )  CYCLE xn
190
191      yn: DO  nya = ny, 2*ny
192!
193!--          Meet conditions for ny
194             IF ( MOD( nya+1, pdims(2) ) /= 0 .OR. &
195                  MOD( nya+1, pdims(1) ) /= 0 )  CYCLE yn
196
197
198         zn: DO  nza = nz, 2*nz
199!
200!--             Meet conditions for nz
201                IF ( ( MOD( nza, pdims(1) ) /= 0  .AND.  pdims(1) /= 1  .AND. &
202                       pdims(2) /= 1 )  .OR.                                  &
203                     ( MOD( nza, pdims(2) ) /= 0  .AND.  dt_dosp /= 9999999.9 &
204                     ) )  THEN
205                   CYCLE zn
206                ELSE
207                   found = .TRUE.
208                   EXIT xn
209                ENDIF
210
211             ENDDO zn
212
213          ENDDO yn
214
215       ENDDO xn
216
217       IF ( .NOT. found )  THEN
218          IF ( myid == 0 )  THEN
219             PRINT*,'+++ init_pegrid: no matching grid for transpositions found'
220          ENDIF
221          CALL local_stop
222       ENDIF
223
224    ENDIF
225
226!
227!-- Calculate array bounds in x-direction for every PE.
228!-- The last PE along x may get less grid points than the others
229    ALLOCATE( nxlf(0:pdims(1)-1), nxrf(0:pdims(1)-1), nynf(0:pdims(2)-1), &
230              nysf(0:pdims(2)-1), nnx_pe(0:pdims(1)-1), nny_pe(0:pdims(2)-1) )
231
232    IF ( MOD( nxa+1 , pdims(1) ) /= 0 )  THEN
233       IF ( myid == 0 )  THEN
234          PRINT*,'+++ x-direction:  gridpoint number (',nx+1,') is not an'
235          PRINT*,'                  integral divisor of the number of proces', &
236                                   &'sors (', pdims(1),')'
237       ENDIF
238       CALL local_stop
239    ELSE
240       nnx  = ( nxa + 1 ) / pdims(1)
241       IF ( nnx*pdims(1) - ( nx + 1) > nnx )  THEN
242          IF ( myid == 0 )  THEN
243             PRINT*,'+++ x-direction: nx does not match the requirements ', &
244                         'given by the number of PEs'
245             PRINT*,'                 used'
246             PRINT*,'    please use nx = ', nx - ( pdims(1) - ( nnx*pdims(1) &
247                         - ( nx + 1 ) ) ), ' instead of nx =', nx
248          ENDIF
249          CALL local_stop
250       ENDIF
251    ENDIF   
252
253!
254!-- Left and right array bounds, number of gridpoints
255    DO  i = 0, pdims(1)-1
256       nxlf(i)   = i * nnx
257       nxrf(i)   = ( i + 1 ) * nnx - 1
258       nnx_pe(i) = MIN( nx, nxrf(i) ) - nxlf(i) + 1
259    ENDDO
260
261!
262!-- Calculate array bounds in y-direction for every PE.
263    IF ( MOD( nya+1 , pdims(2) ) /= 0 )  THEN
264       IF ( myid == 0 )  THEN
265          PRINT*,'+++ y-direction:  gridpoint number (',ny+1,') is not an'
266          PRINT*,'                  integral divisor of the number of proces', &
267                                   &'sors (', pdims(2),')'
268       ENDIF
269       CALL local_stop
270    ELSE
271       nny  = ( nya + 1 ) / pdims(2)
272       IF ( nny*pdims(2) - ( ny + 1) > nny )  THEN
273          IF ( myid == 0 )  THEN
274             PRINT*,'+++ x-direction: nx does not match the requirements ', &
275                         'given by the number of PEs'
276             PRINT*,'                 used'
277             PRINT*,'    please use nx = ', nx - ( pdims(1) - ( nnx*pdims(1) &
278                         - ( nx + 1 ) ) ), ' instead of nx =', nx
279          ENDIF
280          CALL local_stop
281       ENDIF
282    ENDIF   
283
284!
285!-- South and north array bounds
286    DO  j = 0, pdims(2)-1
287       nysf(j)   = j * nny
288       nynf(j)   = ( j + 1 ) * nny - 1
289       nny_pe(j) = MIN( ny, nynf(j) ) - nysf(j) + 1
290    ENDDO
291
292!
293!-- Local array bounds of the respective PEs
294    nxl  = nxlf(pcoord(1))
295    nxra = nxrf(pcoord(1))
296    nxr  = MIN( nx, nxra )
297    nys  = nysf(pcoord(2))
298    nyna = nynf(pcoord(2))
299    nyn  = MIN( ny, nyna )
300    nzb  = 0
301    nzta = nza
302    nzt  = MIN( nz, nzta )
303    nnz  = nza
304
305!
306!-- Calculate array bounds and gridpoint numbers for the transposed arrays
307!-- (needed in the pressure solver)
308!-- For the transposed arrays, cyclic boundaries as well as top and bottom
309!-- boundaries are omitted, because they are obstructive to the transposition
310
311!
312!-- 1. transposition  z --> x
313!-- This transposition is not neccessary in case of a 1d-decomposition along x,
314!-- except that the uptream-spline method is switched on
315    IF ( pdims(2) /= 1  .OR.  momentum_advec == 'ups-scheme'  .OR. &
316         scalar_advec == 'ups-scheme' )  THEN
317
318       IF ( pdims(2) == 1  .AND. ( momentum_advec == 'ups-scheme'  .OR. &
319            scalar_advec == 'ups-scheme' ) )  THEN
320          IF ( myid == 0 )  THEN
321             PRINT*,'+++ WARNING: init_pegrid: 1d-decomposition along x ', &
322                                &'chosen but nz restrictions may occur'
323             PRINT*,'             since ups-scheme is activated'
324          ENDIF
325       ENDIF
326       nys_x  = nys
327       nyn_xa = nyna
328       nyn_x  = nyn
329       nny_x  = nny
330       IF ( MOD( nza , pdims(1) ) /= 0 )  THEN
331          IF ( myid == 0 )  THEN
332             PRINT*,'+++ transposition z --> x:'
333             PRINT*,'    nz=',nz,' is not an integral divisior of pdims(1)=', &
334                    &pdims(1)
335          ENDIF
336          CALL local_stop
337       ENDIF
338       nnz_x  = nza / pdims(1)
339       nzb_x  = 1 + myidx * nnz_x
340       nzt_xa = ( myidx + 1 ) * nnz_x
341       nzt_x  = MIN( nzt, nzt_xa )
342
343       sendrecvcount_zx = nnx * nny * nnz_x
344
345    ELSE
346!
347!---   Setting of dummy values because otherwise variables are undefined in
348!---   the next step  x --> y
349!---   WARNING: This case has still to be clarified!!!!!!!!!!!!
350       nnz_x  = 1
351       nzb_x  = 1
352       nzt_xa = 1
353       nzt_x  = 1
354       nny_x  = nny
355
356    ENDIF
357
358!
359!-- 2. transposition  x --> y
360    nnz_y  = nnz_x
361    nzb_y  = nzb_x
362    nzt_ya = nzt_xa
363    nzt_y  = nzt_x
364    IF ( MOD( nxa+1 , pdims(2) ) /= 0 )  THEN
365       IF ( myid == 0 )  THEN
366          PRINT*,'+++ transposition x --> y:'
367          PRINT*,'    nx+1=',nx+1,' is not an integral divisor of ',&
368                 &'pdims(2)=',pdims(2)
369       ENDIF
370       CALL local_stop
371    ENDIF
372    nnx_y = (nxa+1) / pdims(2)
373    nxl_y = myidy * nnx_y
374    nxr_ya = ( myidy + 1 ) * nnx_y - 1
375    nxr_y  = MIN( nx, nxr_ya )
376
377    sendrecvcount_xy = nnx_y * nny_x * nnz_y
378
379!
380!-- 3. transposition  y --> z  (ELSE:  x --> y  in case of 1D-decomposition
381!-- along x)
382    IF ( pdims(2) /= 1  .OR.  momentum_advec == 'ups-scheme'  .OR. &
383         scalar_advec == 'ups-scheme' )  THEN
384!
385!--    y --> z
386!--    This transposition is not neccessary in case of a 1d-decomposition
387!--    along x, except that the uptream-spline method is switched on
388       nnx_z  = nnx_y
389       nxl_z  = nxl_y
390       nxr_za = nxr_ya
391       nxr_z  = nxr_y
392       IF ( MOD( nya+1 , pdims(1) ) /= 0 )  THEN
393          IF ( myid == 0 )  THEN
394             PRINT*,'+++ Transposition y --> z:'
395             PRINT*,'    ny+1=',ny+1,' is not an integral divisor of ',&
396                    &'pdims(1)=',pdims(1)
397          ENDIF
398          CALL local_stop
399       ENDIF
400       nny_z  = (nya+1) / pdims(1)
401       nys_z  = myidx * nny_z
402       nyn_za = ( myidx + 1 ) * nny_z - 1
403       nyn_z  = MIN( ny, nyn_za )
404
405       sendrecvcount_yz = nnx_y * nny_z * nnz_y
406
407    ELSE
408!
409!--    x --> y. This condition must be fulfilled for a 1D-decomposition along x
410       IF ( MOD( nya+1 , pdims(1) ) /= 0 )  THEN
411          IF ( myid == 0 )  THEN
412             PRINT*,'+++ Transposition x --> y:'
413             PRINT*,'    ny+1=',ny+1,' is not an integral divisor of ',&
414                    &'pdims(1)=',pdims(1)
415          ENDIF
416          CALL local_stop
417       ENDIF
418
419    ENDIF
420
421!
422!-- Indices for direct transpositions z --> y (used for calculating spectra)
423    IF ( dt_dosp /= 9999999.9 )  THEN
424       IF ( MOD( nza, pdims(2) ) /= 0 )  THEN
425          IF ( myid == 0 )  THEN
426             PRINT*,'+++ Direct transposition z --> y (needed for spectra):'
427             PRINT*,'    nz=',nz,' is not an integral divisor of ',&
428                    &'pdims(2)=',pdims(2)
429          ENDIF
430          CALL local_stop
431       ELSE
432          nxl_yd  = nxl
433          nxr_yda = nxra
434          nxr_yd  = nxr
435          nzb_yd  = 1 + myidy * ( nza / pdims(2) )
436          nzt_yda = ( myidy + 1 ) * ( nza / pdims(2) )
437          nzt_yd  = MIN( nzt, nzt_yda )
438
439          sendrecvcount_zyd = nnx * nny * ( nza / pdims(2) )
440       ENDIF
441    ENDIF
442
443!
444!-- Indices for direct transpositions y --> x (they are only possible in case
445!-- of a 1d-decomposition along x)
446    IF ( pdims(2) == 1 )  THEN
447       nny_x  = nny / pdims(1)
448       nys_x  = myid * nny_x
449       nyn_xa = ( myid + 1 ) * nny_x - 1
450       nyn_x  = MIN( ny, nyn_xa )
451       nzb_x  = 1
452       nzt_xa = nza
453       nzt_x  = nz
454       sendrecvcount_xy = nnx * nny_x * nza
455    ENDIF
456
457!
458!-- Indices for direct transpositions x --> y (they are only possible in case
459!-- of a 1d-decomposition along y)
460    IF ( pdims(1) == 1 )  THEN
461       nnx_y  = nnx / pdims(2)
462       nxl_y  = myid * nnx_y
463       nxr_ya = ( myid + 1 ) * nnx_y - 1
464       nxr_y  = MIN( nx, nxr_ya )
465       nzb_y  = 1
466       nzt_ya = nza
467       nzt_y  = nz
468       sendrecvcount_xy = nnx_y * nny * nza
469    ENDIF
470
471!
472!-- Arrays for storing the array bounds are needed any more
473    DEALLOCATE( nxlf , nxrf , nynf , nysf )
474
475!
476!-- Collect index bounds from other PEs (to be written to restart file later)
477    ALLOCATE( hor_index_bounds(4,0:numprocs-1) )
478
479    IF ( myid == 0 )  THEN
480
481       hor_index_bounds(1,0) = nxl
482       hor_index_bounds(2,0) = nxr
483       hor_index_bounds(3,0) = nys
484       hor_index_bounds(4,0) = nyn
485
486!
487!--    Receive data from all other PEs
488       DO  i = 1, numprocs-1
489          CALL MPI_RECV( ibuf, 4, MPI_INTEGER, i, MPI_ANY_TAG, comm2d, status, &
490                         ierr )
491          hor_index_bounds(:,i) = ibuf(1:4)
492       ENDDO
493
494    ELSE
495!
496!--    Send index bounds to PE0
497       ibuf(1) = nxl
498       ibuf(2) = nxr
499       ibuf(3) = nys
500       ibuf(4) = nyn
501       CALL MPI_SEND( ibuf, 4, MPI_INTEGER, 0, myid, comm2d, ierr )
502
503    ENDIF
504
505#if defined( __print )
506!
507!-- Control output
508    IF ( myid == 0 )  THEN
509       PRINT*, '*** processor topology ***'
510       PRINT*, ' '
511       PRINT*, 'myid   pcoord    left right  south north  idx idy   nxl: nxr',&
512               &'   nys: nyn'
513       PRINT*, '------------------------------------------------------------',&
514               &'-----------'
515       WRITE (*,1000)  0, pcoord(1), pcoord(2), pleft, pright, psouth, pnorth, &
516                       myidx, myidy, nxl, nxr, nys, nyn
5171000   FORMAT (I4,2X,'(',I3,',',I3,')',3X,I4,2X,I4,3X,I4,2X,I4,2X,I3,1X,I3, &
518               2(2X,I4,':',I4))
519
520!
521!--    Receive data from the other PEs
522       DO  i = 1,numprocs-1
523          CALL MPI_RECV( ibuf, 12, MPI_INTEGER, i, MPI_ANY_TAG, comm2d, status, &
524                         ierr )
525          WRITE (*,1000)  i, ( ibuf(j) , j = 1,12 )
526       ENDDO
527    ELSE
528
529!
530!--    Send data to PE0
531       ibuf(1) = pcoord(1); ibuf(2) = pcoord(2); ibuf(3) = pleft
532       ibuf(4) = pright; ibuf(5) = psouth; ibuf(6) = pnorth; ibuf(7) = myidx
533       ibuf(8) = myidy; ibuf(9) = nxl; ibuf(10) = nxr; ibuf(11) = nys
534       ibuf(12) = nyn
535       CALL MPI_SEND( ibuf, 12, MPI_INTEGER, 0, myid, comm2d, ierr )       
536    ENDIF
537#endif
538
539#if defined( __mpi2 )
540!
541!-- In case of coupled runs, get the port name on PE0 of the atmosphere model
542!-- and pass it to PE0 of the ocean model
543    IF ( myid == 0 )  THEN
544
545       IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
546
547          CALL MPI_OPEN_PORT( MPI_INFO_NULL, port_name, ierr )
548!
549!--       TEST OUTPUT (TO BE REMOVED)
550          WRITE(9,*)  TRIM( coupling_mode ),  &
551               ', ierr after MPI_OPEN_PORT: ', ierr
552          CALL LOCAL_FLUSH( 9 )
553
554          CALL MPI_PUBLISH_NAME( 'palm_coupler', MPI_INFO_NULL, port_name, &
555                                 ierr )
556!
557!--       TEST OUTPUT (TO BE REMOVED)
558          WRITE(9,*)  TRIM( coupling_mode ),  &
559               ', ierr after MPI_PUBLISH_NAME: ', ierr
560          CALL LOCAL_FLUSH( 9 )
561
562!
563!--       Write a flag file for the ocean model and the other atmosphere
564!--       processes.
565!--       There seems to be a bug in MPICH2 which causes hanging processes
566!--       in case that execution of LOOKUP_NAME is continued too early
567!--       (i.e. before the port has been created)
568          OPEN( 90, FILE='COUPLING_PORT_OPENED', FORM='FORMATTED' )
569          WRITE ( 90, '(''TRUE'')' )
570          CLOSE ( 90 )
571
572       ELSEIF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
573
574!
575!--       Continue only if the atmosphere model has created the port.
576!--       There seems to be a bug in MPICH2 which causes hanging processes
577!--       in case that execution of LOOKUP_NAME is continued too early
578!--       (i.e. before the port has been created)
579          INQUIRE( FILE='COUPLING_PORT_OPENED', EXIST=found )
580          DO WHILE ( .NOT. found )
581             INQUIRE( FILE='COUPLING_PORT_OPENED', EXIST=found )
582          ENDDO
583
584          CALL MPI_LOOKUP_NAME( 'palm_coupler', MPI_INFO_NULL, port_name, ierr )
585!
586!--       TEST OUTPUT (TO BE REMOVED)
587          WRITE(9,*)  TRIM( coupling_mode ),  &
588               ', ierr after MPI_LOOKUP_NAME: ', ierr
589          CALL LOCAL_FLUSH( 9 )
590
591
592       ENDIF
593
594    ENDIF
595
596!
597!-- In case of coupled runs, establish the connection between the atmosphere
598!-- and the ocean model and define the intercommunicator (comm_inter)
599    CALL MPI_BARRIER( comm2d, ierr )
600    IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
601
602       print*, '... before COMM_ACCEPT'
603       CALL MPI_COMM_ACCEPT( port_name, MPI_INFO_NULL, 0, MPI_COMM_WORLD, &
604                             comm_inter, ierr )
605       print*, '--- ierr = ', ierr
606       print*, '--- comm_inter atmosphere = ', comm_inter
607
608       coupling_mode_remote = 'ocean_to_atmosphere'
609
610    ELSEIF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
611
612       IF ( myid == 0 )  PRINT*, '*** read: ', port_name, '  ierr = ', ierr
613       print*, '... before COMM_CONNECT'
614       CALL MPI_COMM_CONNECT( port_name, MPI_INFO_NULL, 0, MPI_COMM_WORLD, &
615                              comm_inter, ierr )
616       print*, '--- ierr = ', ierr
617       print*, '--- comm_inter ocean      = ', comm_inter
618
619       coupling_mode_remote = 'atmosphere_to_ocean'
620
621    ENDIF
622
623!
624!-- In case of coupled runs, create a new MPI derived datatype for the
625!-- exchange of surface (xy) data .
626!-- Gridpoint number for the exchange of ghost points (xy-plane)
627    ngp_xy  = ( nxr - nxl + 3 ) * ( nyn - nys + 3 )
628
629!
630!-- Define a new MPI derived datatype for the exchange of ghost points in
631!-- y-direction for 2D-arrays (line)
632    CALL MPI_TYPE_VECTOR( ngp_xy, 1, nzt-nzb+2, MPI_REAL, type_xy, ierr )
633    CALL MPI_TYPE_COMMIT( type_xy, ierr )
634#endif
635
636#else
637
638!
639!-- Array bounds when running on a single PE (respectively a non-parallel
640!-- machine)
641    nxl  = 0
642    nxr  = nx
643    nxra = nx
644    nnx  = nxr - nxl + 1
645    nys  = 0
646    nyn  = ny
647    nyna = ny
648    nny  = nyn - nys + 1
649    nzb  = 0
650    nzt  = nz
651    nzta = nz
652    nnz  = nz
653
654    ALLOCATE( hor_index_bounds(4,0:0) )
655    hor_index_bounds(1,0) = nxl
656    hor_index_bounds(2,0) = nxr
657    hor_index_bounds(3,0) = nys
658    hor_index_bounds(4,0) = nyn
659
660!
661!-- Array bounds for the pressure solver (in the parallel code, these bounds
662!-- are the ones for the transposed arrays)
663    nys_x  = nys
664    nyn_x  = nyn
665    nyn_xa = nyn
666    nzb_x  = nzb + 1
667    nzt_x  = nzt
668    nzt_xa = nzt
669
670    nxl_y  = nxl
671    nxr_y  = nxr
672    nxr_ya = nxr
673    nzb_y  = nzb + 1
674    nzt_y  = nzt
675    nzt_ya = nzt
676
677    nxl_z  = nxl
678    nxr_z  = nxr
679    nxr_za = nxr
680    nys_z  = nys
681    nyn_z  = nyn
682    nyn_za = nyn
683
684#endif
685
686!
687!-- Calculate number of grid levels necessary for the multigrid poisson solver
688!-- as well as the gridpoint indices on each level
689    IF ( psolver == 'multigrid' )  THEN
690
691!
692!--    First calculate number of possible grid levels for the subdomains
693       mg_levels_x = 1
694       mg_levels_y = 1
695       mg_levels_z = 1
696
697       i = nnx
698       DO WHILE ( MOD( i, 2 ) == 0  .AND.  i /= 2 )
699          i = i / 2
700          mg_levels_x = mg_levels_x + 1
701       ENDDO
702
703       j = nny
704       DO WHILE ( MOD( j, 2 ) == 0  .AND.  j /= 2 )
705          j = j / 2
706          mg_levels_y = mg_levels_y + 1
707       ENDDO
708
709       k = nz    ! do not use nnz because it might be > nz due to transposition
710                 ! requirements
711       DO WHILE ( MOD( k, 2 ) == 0  .AND.  k /= 2 )
712          k = k / 2
713          mg_levels_z = mg_levels_z + 1
714       ENDDO
715
716       maximum_grid_level = MIN( mg_levels_x, mg_levels_y, mg_levels_z )
717
718!
719!--    Find out, if the total domain allows more levels. These additional
720!--    levels are processed on PE0 only.
721       IF ( numprocs > 1 )  THEN
722          IF ( mg_levels_z > MIN( mg_levels_x, mg_levels_y ) )  THEN
723             mg_switch_to_pe0_level_l = maximum_grid_level
724
725             mg_levels_x = 1
726             mg_levels_y = 1
727
728             i = nx+1
729             DO WHILE ( MOD( i, 2 ) == 0  .AND.  i /= 2 )
730                i = i / 2
731                mg_levels_x = mg_levels_x + 1
732             ENDDO
733
734             j = ny+1
735             DO WHILE ( MOD( j, 2 ) == 0  .AND.  j /= 2 )
736                j = j / 2
737                mg_levels_y = mg_levels_y + 1
738             ENDDO
739
740             maximum_grid_level_l = MIN( mg_levels_x, mg_levels_y, mg_levels_z )
741
742             IF ( maximum_grid_level_l > mg_switch_to_pe0_level_l )  THEN
743                mg_switch_to_pe0_level_l = maximum_grid_level_l - &
744                                           mg_switch_to_pe0_level_l + 1
745             ELSE
746                mg_switch_to_pe0_level_l = 0
747             ENDIF
748          ELSE
749             mg_switch_to_pe0_level_l = 0
750             maximum_grid_level_l = maximum_grid_level
751          ENDIF
752
753!
754!--       Use switch level calculated above only if it is not pre-defined
755!--       by user
756          IF ( mg_switch_to_pe0_level == 0 )  THEN
757
758             IF ( mg_switch_to_pe0_level_l /= 0 )  THEN
759                mg_switch_to_pe0_level = mg_switch_to_pe0_level_l
760                maximum_grid_level     = maximum_grid_level_l
761             ENDIF
762
763          ELSE
764!
765!--          Check pre-defined value and reset to default, if neccessary
766             IF ( mg_switch_to_pe0_level < mg_switch_to_pe0_level_l  .OR.  &
767                  mg_switch_to_pe0_level >= maximum_grid_level_l )  THEN
768                IF ( myid == 0 )  THEN
769                   PRINT*, '+++ WARNING init_pegrid: mg_switch_to_pe0_level ', &
770                               'out of range and reset to default (=0)'
771                ENDIF
772                mg_switch_to_pe0_level = 0
773             ELSE
774!
775!--             Use the largest number of possible levels anyway and recalculate
776!--             the switch level to this largest number of possible values
777                maximum_grid_level = maximum_grid_level_l
778
779             ENDIF
780          ENDIF
781
782       ENDIF
783
784       ALLOCATE( grid_level_count(maximum_grid_level),                   &
785                 nxl_mg(maximum_grid_level), nxr_mg(maximum_grid_level), &
786                 nyn_mg(maximum_grid_level), nys_mg(maximum_grid_level), &
787                 nzt_mg(maximum_grid_level) )
788
789       grid_level_count = 0
790       nxl_l = nxl; nxr_l = nxr; nys_l = nys; nyn_l = nyn; nzt_l = nzt
791
792       DO  i = maximum_grid_level, 1 , -1
793
794          IF ( i == mg_switch_to_pe0_level )  THEN
795#if defined( __parallel )
796!
797!--          Save the grid size of the subdomain at the switch level, because
798!--          it is needed in poismg.
799!--          Array bounds of the local subdomain grids are gathered on PE0
800             ind(1) = nxl_l; ind(2) = nxr_l
801             ind(3) = nys_l; ind(4) = nyn_l
802             ind(5) = nzt_l
803             ALLOCATE( ind_all(5*numprocs), mg_loc_ind(5,0:numprocs-1) )
804             CALL MPI_ALLGATHER( ind, 5, MPI_INTEGER, ind_all, 5, &
805                                 MPI_INTEGER, comm2d, ierr )
806             DO  j = 0, numprocs-1
807                DO  k = 1, 5
808                   mg_loc_ind(k,j) = ind_all(k+j*5)
809                ENDDO
810             ENDDO
811             DEALLOCATE( ind_all )
812!
813!--          Calculate the grid size of the total domain gathered on PE0
814             nxr_l = ( nxr_l-nxl_l+1 ) * pdims(1) - 1
815             nxl_l = 0
816             nyn_l = ( nyn_l-nys_l+1 ) * pdims(2) - 1
817             nys_l = 0
818!
819!--          The size of this gathered array must not be larger than the
820!--          array tend, which is used in the multigrid scheme as a temporary
821!--          array
822             subdomain_size = ( nxr - nxl + 3 )     * ( nyn - nys + 3 )     * &
823                              ( nzt - nzb + 2 )
824             gathered_size  = ( nxr_l - nxl_l + 3 ) * ( nyn_l - nys_l + 3 ) * &
825                              ( nzt_l - nzb + 2 )
826
827             IF ( gathered_size > subdomain_size )  THEN
828                IF ( myid == 0 )  THEN
829                   PRINT*, '+++ init_pegrid: not enough memory for storing ', &
830                               'gathered multigrid data on PE0'
831                ENDIF
832                CALL local_stop
833             ENDIF
834#else
835             PRINT*, '+++ init_pegrid: multigrid gather/scatter impossible ', &
836                          'in non parallel mode'
837             CALL local_stop
838#endif
839          ENDIF
840
841          nxl_mg(i) = nxl_l
842          nxr_mg(i) = nxr_l
843          nys_mg(i) = nys_l
844          nyn_mg(i) = nyn_l
845          nzt_mg(i) = nzt_l
846
847          nxl_l = nxl_l / 2 
848          nxr_l = nxr_l / 2
849          nys_l = nys_l / 2 
850          nyn_l = nyn_l / 2 
851          nzt_l = nzt_l / 2 
852       ENDDO
853
854    ELSE
855
856       maximum_grid_level = 1
857
858    ENDIF
859
860    grid_level = maximum_grid_level
861
862#if defined( __parallel )
863!
864!-- Gridpoint number for the exchange of ghost points (y-line for 2D-arrays)
865    ngp_y  = nyn - nys + 1
866
867!
868!-- Define a new MPI derived datatype for the exchange of ghost points in
869!-- y-direction for 2D-arrays (line)
870    CALL MPI_TYPE_VECTOR( nxr-nxl+3, 1, ngp_y+2, MPI_REAL, type_x, ierr )
871    CALL MPI_TYPE_COMMIT( type_x, ierr )
872    CALL MPI_TYPE_VECTOR( nxr-nxl+3, 1, ngp_y+2, MPI_INTEGER, type_x_int, ierr )
873    CALL MPI_TYPE_COMMIT( type_x_int, ierr )
874
875!
876!-- Calculate gridpoint numbers for the exchange of ghost points along x
877!-- (yz-plane for 3D-arrays) and define MPI derived data type(s) for the
878!-- exchange of ghost points in y-direction (xz-plane).
879!-- Do these calculations for the model grid and (if necessary) also
880!-- for the coarser grid levels used in the multigrid method
881    ALLOCATE ( ngp_yz(maximum_grid_level), type_xz(maximum_grid_level) )
882
883    nxl_l = nxl; nxr_l = nxr; nys_l = nys; nyn_l = nyn; nzb_l = nzb; nzt_l = nzt
884         
885    DO i = maximum_grid_level, 1 , -1
886       ngp_yz(i) = (nzt_l - nzb_l + 2) * (nyn_l - nys_l + 3)
887
888       CALL MPI_TYPE_VECTOR( nxr_l-nxl_l+3, nzt_l-nzb_l+2, ngp_yz(i), &
889                             MPI_REAL, type_xz(i), ierr )
890       CALL MPI_TYPE_COMMIT( type_xz(i), ierr )
891
892       nxl_l = nxl_l / 2 
893       nxr_l = nxr_l / 2
894       nys_l = nys_l / 2 
895       nyn_l = nyn_l / 2 
896       nzt_l = nzt_l / 2 
897    ENDDO
898#endif
899
900#if defined( __parallel )
901!
902!-- Setting of flags for inflow/outflow conditions in case of non-cyclic
903!-- horizontal boundary conditions.
904    IF ( pleft == MPI_PROC_NULL )  THEN
905       IF ( bc_lr == 'dirichlet/radiation' )  THEN
906          inflow_l  = .TRUE.
907       ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
908          outflow_l = .TRUE.
909       ENDIF
910    ENDIF
911
912    IF ( pright == MPI_PROC_NULL )  THEN
913       IF ( bc_lr == 'dirichlet/radiation' )  THEN
914          outflow_r = .TRUE.
915       ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
916          inflow_r  = .TRUE.
917       ENDIF
918    ENDIF
919
920    IF ( psouth == MPI_PROC_NULL )  THEN
921       IF ( bc_ns == 'dirichlet/radiation' )  THEN
922          outflow_s = .TRUE.
923       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
924          inflow_s  = .TRUE.
925       ENDIF
926    ENDIF
927
928    IF ( pnorth == MPI_PROC_NULL )  THEN
929       IF ( bc_ns == 'dirichlet/radiation' )  THEN
930          inflow_n  = .TRUE.
931       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
932          outflow_n = .TRUE.
933       ENDIF
934    ENDIF
935
936!
937!-- Broadcast the id of the inflow PE
938    IF ( inflow_l )  THEN
939       id_inflow_l = myidx
940    ELSE
941       id_inflow_l = 0
942    ENDIF
943    CALL MPI_ALLREDUCE( id_inflow_l, id_inflow, 1, MPI_INTEGER, MPI_SUM, &
944                        comm1dx, ierr )
945
946!
947!-- Broadcast the id of the recycling plane
948!-- WARNING: needs to be adjusted in case of inflows other than from left side!
949    IF ( ( recycling_width / dx ) >= nxl  .AND.  ( recycling_width / dx ) <= nxr ) &
950    THEN
951       id_recycling_l = myidx
952    ELSE
953       id_recycling_l = 0
954    ENDIF
955    CALL MPI_ALLREDUCE( id_recycling_l, id_recycling, 1, MPI_INTEGER, MPI_SUM, &
956                        comm1dx, ierr )
957
958#else
959    IF ( bc_lr == 'dirichlet/radiation' )  THEN
960       inflow_l  = .TRUE.
961       outflow_r = .TRUE.
962    ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
963       outflow_l = .TRUE.
964       inflow_r  = .TRUE.
965    ENDIF
966
967    IF ( bc_ns == 'dirichlet/radiation' )  THEN
968       inflow_n  = .TRUE.
969       outflow_s = .TRUE.
970    ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
971       outflow_n = .TRUE.
972       inflow_s  = .TRUE.
973    ENDIF
974#endif
975!
976!-- At the outflow, u or v, respectively, have to be calculated for one more
977!-- grid point.
978    IF ( outflow_l )  THEN
979       nxlu = nxl + 1
980    ELSE
981       nxlu = nxl
982    ENDIF
983    IF ( outflow_s )  THEN
984       nysv = nys + 1
985    ELSE
986       nysv = nys
987    ENDIF
988
989    IF ( psolver == 'poisfft_hybrid' )  THEN
990       CALL poisfft_hybrid_ini
991    ELSEIF ( psolver == 'poisfft' )  THEN
992       CALL poisfft_init
993    ENDIF
994
995!
996!-- Allocate wall flag arrays used in the multigrid solver
997    IF ( psolver == 'multigrid' )  THEN
998
999       DO  i = maximum_grid_level, 1, -1
1000
1001           SELECT CASE ( i )
1002
1003              CASE ( 1 )
1004                 ALLOCATE( wall_flags_1(nzb:nzt_mg(i)+1,         &
1005                                        nys_mg(i)-1:nyn_mg(i)+1, &
1006                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1007
1008              CASE ( 2 )
1009                 ALLOCATE( wall_flags_2(nzb:nzt_mg(i)+1,         &
1010                                        nys_mg(i)-1:nyn_mg(i)+1, &
1011                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1012
1013              CASE ( 3 )
1014                 ALLOCATE( wall_flags_3(nzb:nzt_mg(i)+1,         &
1015                                        nys_mg(i)-1:nyn_mg(i)+1, &
1016                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1017
1018              CASE ( 4 )
1019                 ALLOCATE( wall_flags_4(nzb:nzt_mg(i)+1,         &
1020                                        nys_mg(i)-1:nyn_mg(i)+1, &
1021                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1022
1023              CASE ( 5 )
1024                 ALLOCATE( wall_flags_5(nzb:nzt_mg(i)+1,         &
1025                                        nys_mg(i)-1:nyn_mg(i)+1, &
1026                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1027
1028              CASE ( 6 )
1029                 ALLOCATE( wall_flags_6(nzb:nzt_mg(i)+1,         &
1030                                        nys_mg(i)-1:nyn_mg(i)+1, &
1031                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1032
1033              CASE ( 7 )
1034                 ALLOCATE( wall_flags_7(nzb:nzt_mg(i)+1,         &
1035                                        nys_mg(i)-1:nyn_mg(i)+1, &
1036                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1037
1038              CASE ( 8 )
1039                 ALLOCATE( wall_flags_8(nzb:nzt_mg(i)+1,         &
1040                                        nys_mg(i)-1:nyn_mg(i)+1, &
1041                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1042
1043              CASE ( 9 )
1044                 ALLOCATE( wall_flags_9(nzb:nzt_mg(i)+1,         &
1045                                        nys_mg(i)-1:nyn_mg(i)+1, &
1046                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1047
1048              CASE ( 10 )
1049                 ALLOCATE( wall_flags_10(nzb:nzt_mg(i)+1,        &
1050                                        nys_mg(i)-1:nyn_mg(i)+1, &
1051                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1052
1053              CASE DEFAULT
1054                 IF ( myid == 0 )  PRINT*, '+++ init_pegrid: more than 10 ', &
1055                                           ' multigrid levels'
1056                 CALL local_stop
1057
1058          END SELECT
1059
1060       ENDDO
1061
1062    ENDIF
1063
1064 END SUBROUTINE init_pegrid
Note: See TracBrowser for help on using the repository browser.