Home

Context Navigation

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

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |