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

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

bugfixes for turbulent inflow in init_pegrid, inflow_turbulence, and init_3d_model

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