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

Last change on this file since 647 was 647, checked in by raasch, 11 years ago

last commit documented

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