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

Last change on this file since 232 was 226, checked in by raasch, 16 years ago

preparations for the next release

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