Home

Context Navigation

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

Last change on this file since 1374 was 1354, checked in by heinze, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 40.6 KB

Rev	Line
[1]	1	SUBROUTINE init_pegrid
[1036]	2
	3	!--------------------------------------------------------------------------------!
	4	! This file is part of PALM.
	5	!
	6	! PALM is free software: you can redistribute it and/or modify it under the terms
	7	! of the GNU General Public License as published by the Free Software Foundation,
	8	! either version 3 of the License, or (at your option) any later version.
	9	!
	10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
	11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
	12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
	13	!
	14	! You should have received a copy of the GNU General Public License along with
	15	! PALM. If not, see <http://www.gnu.org/licenses/>.
	16	!
[1310]	17	! Copyright 1997-2014 Leibniz Universitaet Hannover
[1036]	18	!--------------------------------------------------------------------------------!
	19	!
[254]	20	! Current revisions:
[1322]	21	! ------------------
[1354]	22	!
	23	!
[1321]	24	! Former revisions:
	25	! -----------------
	26	! $Id: init_pegrid.f90 1354 2014-04-08 15:22:57Z raasch $
	27	!
[1354]	28	! 1353 2014-04-08 15:21:23Z heinze
	29	! REAL constants provided with KIND-attribute
	30	!
[1323]	31	! 1322 2014-03-20 16:38:49Z raasch
	32	! REAL functions provided with KIND-attribute
	33	!
[1321]	34	! 1320 2014-03-20 08:40:49Z raasch
[1320]	35	! ONLY-attribute added to USE-statements,
	36	! kind-parameters added to all INTEGER and REAL declaration statements,
	37	! kinds are defined in new module kinds,
	38	! revision history before 2012 removed,
	39	! comment fields (!:) to be used for variable explanations added to
	40	! all variable declaration statements
[760]	41	!
[1305]	42	! 1304 2014-03-12 10:29:42Z raasch
	43	! bugfix: single core MPI runs missed some settings of transpose indices
	44	!
[1213]	45	! 1212 2013-08-15 08:46:27Z raasch
	46	! error message for poisfft_hybrid removed
	47	!
[1160]	48	! 1159 2013-05-21 11:58:22Z fricke
	49	! dirichlet/neumann and neumann/dirichlet removed
	50	!
[1140]	51	! 1139 2013-04-18 07:25:03Z raasch
	52	! bugfix for calculating the id of the PE carrying the recycling plane
	53	!
[1112]	54	! 1111 2013-03-08 23:54:10Z raasch
	55	! initialization of poisfft moved to module poisfft
	56	!
[1093]	57	! 1092 2013-02-02 11:24:22Z raasch
	58	! unused variables removed
	59	!
[1057]	60	! 1056 2012-11-16 15:28:04Z raasch
	61	! Indices for arrays n.._mg start from zero due to definition of arrays f2 and
	62	! p2 as automatic arrays in recursive subroutine next_mg_level
	63	!
[1042]	64	! 1041 2012-11-06 02:36:29Z raasch
	65	! a 2d virtual processor topology is used by default for all machines
	66	!
[1037]	67	! 1036 2012-10-22 13:43:42Z raasch
	68	! code put under GPL (PALM 3.9)
	69	!
[1004]	70	! 1003 2012-09-14 14:35:53Z raasch
	71	! subdomains must have identical size (grid matching = "match" removed)
	72	!
[1002]	73	! 1001 2012-09-13 14:08:46Z raasch
	74	! all actions concerning upstream-spline-method removed
	75	!
[979]	76	! 978 2012-08-09 08:28:32Z fricke
	77	! dirichlet/neumann and neumann/dirichlet added
	78	! nxlu and nysv are also calculated for inflow boundary
	79	!
[810]	80	! 809 2012-01-30 13:32:58Z maronga
	81	! Bugfix: replaced .AND. and .NOT. with && and ! in the preprocessor directives
	82	!
[808]	83	! 807 2012-01-25 11:53:51Z maronga
	84	! New cpp directive "__check" implemented which is used by check_namelist_files
	85	!
[1]	86	! Revision 1.1 1997/07/24 11:15:09 raasch
	87	! Initial revision
	88	!
	89	!
	90	! Description:
	91	! ------------
	92	! Determination of the virtual processor topology (if not prescribed by the
	93	! user)and computation of the grid point number and array bounds of the local
	94	! domains.
	95	!------------------------------------------------------------------------------!
	96
[1320]	97	USE control_parameters, &
	98	ONLY: bc_lr, bc_ns, coupling_mode, coupling_topology, dt_dosp, &
	99	gathered_size, grid_level, grid_level_count, host, inflow_l, &
	100	inflow_n, inflow_r, inflow_s, io_blocks, io_group, &
	101	maximum_grid_level, maximum_parallel_io_streams, message_string,&
	102	mg_switch_to_pe0_level, momentum_advec, psolver, outflow_l, &
	103	outflow_n, outflow_r, outflow_s, recycling_width, scalar_advec, &
	104	subdomain_size
[1]	105
[1320]	106	USE grid_variables, &
	107	ONLY: dx
	108
	109	USE indices, &
	110	ONLY: mg_loc_ind, nbgp, nnx, nny, nnz, nx, nx_a, nx_o, nxl, nxl_mg, &
	111	nxlu, nxr, nxr_mg, ny, ny_a, ny_o, nyn, nyn_mg, nys, nys_mg, &
	112	nysv, nz, nzb, nzt, nzt_mg, wall_flags_1, wall_flags_2, &
	113	wall_flags_3, wall_flags_4, wall_flags_5, wall_flags_6, &
	114	wall_flags_7, wall_flags_8, wall_flags_9, wall_flags_10
[1]	115
[1320]	116	USE kinds
	117
	118	USE pegrid
	119
	120	USE transpose_indices, &
	121	ONLY: nxl_y, nxl_yd, nxl_z, nxr_y, nxr_yd, nxr_z, nyn_x, nyn_z, nys_x,&
	122	nys_z, nzb_x, nzb_y, nzb_yd, nzt_x, nzt_yd, nzt_y
[667]	123
[1]	124	IMPLICIT NONE
	125
[1320]	126	INTEGER(iwp) :: i !:
	127	INTEGER(iwp) :: id_inflow_l !:
	128	INTEGER(iwp) :: id_recycling_l !:
	129	INTEGER(iwp) :: ind(5) !:
	130	INTEGER(iwp) :: j !:
	131	INTEGER(iwp) :: k !:
	132	INTEGER(iwp) :: maximum_grid_level_l !:
	133	INTEGER(iwp) :: mg_levels_x !:
	134	INTEGER(iwp) :: mg_levels_y !:
	135	INTEGER(iwp) :: mg_levels_z !:
	136	INTEGER(iwp) :: mg_switch_to_pe0_level_l !:
	137	INTEGER(iwp) :: nnx_y !:
	138	INTEGER(iwp) :: nnx_z !:
	139	INTEGER(iwp) :: nny_x !:
	140	INTEGER(iwp) :: nny_z !:
	141	INTEGER(iwp) :: nnz_x !:
	142	INTEGER(iwp) :: nnz_y !:
	143	INTEGER(iwp) :: numproc_sqr !:
	144	INTEGER(iwp) :: nxl_l !:
	145	INTEGER(iwp) :: nxr_l !:
	146	INTEGER(iwp) :: nyn_l !:
	147	INTEGER(iwp) :: nys_l !:
	148	INTEGER(iwp) :: nzb_l !:
	149	INTEGER(iwp) :: nzt_l !:
	150	INTEGER(iwp) :: omp_get_num_threads !:
[1]	151
[1320]	152	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: ind_all !:
	153	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nxlf !:
	154	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nxrf !:
	155	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nynf !:
	156	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nysf !:
[1]	157
[1320]	158	INTEGER(iwp), DIMENSION(2) :: pdims_remote !:
[667]	159
[1092]	160	#if defined( __mpi2 )
[1320]	161	LOGICAL :: found !:
[1092]	162	#endif
[1]	163
	164	!
	165	!-- Get the number of OpenMP threads
	166	!$OMP PARALLEL
[82]	167	#if defined( __intel_openmp_bug )
[1]	168	threads_per_task = omp_get_num_threads()
	169	#else
	170	!$ threads_per_task = omp_get_num_threads()
	171	#endif
	172	!$OMP END PARALLEL
	173
	174
	175	#if defined( __parallel )
[667]	176
[1]	177	!
	178	!-- Determine the processor topology or check it, if prescribed by the user
	179	IF ( npex == -1 .AND. npey == -1 ) THEN
	180
	181	!
	182	!-- Automatic determination of the topology
[1322]	183	numproc_sqr = SQRT( REAL( numprocs, KIND=wp ) )
[1041]	184	pdims(1) = MAX( numproc_sqr , 1 )
	185	DO WHILE ( MOD( numprocs , pdims(1) ) /= 0 )
	186	pdims(1) = pdims(1) - 1
	187	ENDDO
	188	pdims(2) = numprocs / pdims(1)
[1]	189
	190	ELSEIF ( npex /= -1 .AND. npey /= -1 ) THEN
	191
	192	!
	193	!-- Prescribed by user. Number of processors on the prescribed topology
	194	!-- must be equal to the number of PEs available to the job
	195	IF ( ( npex * npey ) /= numprocs ) THEN
[274]	196	WRITE( message_string, * ) 'number of PEs of the prescribed ', &
	197	'topology (', npex*npey,') does not match & the number of ', &
	198	'PEs available to the job (', numprocs, ')'
[254]	199	CALL message( 'init_pegrid', 'PA0221', 1, 2, 0, 6, 0 )
[1]	200	ENDIF
	201	pdims(1) = npex
	202	pdims(2) = npey
	203
	204	ELSE
	205	!
	206	!-- If the processor topology is prescribed by the user, the number of
	207	!-- PEs must be given in both directions
[274]	208	message_string = 'if the processor topology is prescribed by the, ' // &
	209	' user& both values of "npex" and "npey" must be given ' // &
	210	'in the &NAMELIST-parameter file'
[254]	211	CALL message( 'init_pegrid', 'PA0222', 1, 2, 0, 6, 0 )
[1]	212
	213	ENDIF
	214
	215	!
[622]	216	!-- For communication speedup, set barriers in front of collective
	217	!-- communications by default on SGI-type systems
	218	IF ( host(3:5) == 'sgi' ) collective_wait = .TRUE.
	219
	220	!
[1]	221	!-- If necessary, set horizontal boundary conditions to non-cyclic
[722]	222	IF ( bc_lr /= 'cyclic' ) cyclic(1) = .FALSE.
	223	IF ( bc_ns /= 'cyclic' ) cyclic(2) = .FALSE.
[1]	224
[807]	225
[809]	226	#if ! defined( __check)
[1]	227	!
	228	!-- Create the virtual processor grid
	229	CALL MPI_CART_CREATE( comm_palm, ndim, pdims, cyclic, reorder, &
	230	comm2d, ierr )
	231	CALL MPI_COMM_RANK( comm2d, myid, ierr )
	232	WRITE (myid_char,'(''_'',I4.4)') myid
	233
	234	CALL MPI_CART_COORDS( comm2d, myid, ndim, pcoord, ierr )
	235	CALL MPI_CART_SHIFT( comm2d, 0, 1, pleft, pright, ierr )
	236	CALL MPI_CART_SHIFT( comm2d, 1, 1, psouth, pnorth, ierr )
	237
	238	!
	239	!-- Determine sub-topologies for transpositions
	240	!-- Transposition from z to x:
	241	remain_dims(1) = .TRUE.
	242	remain_dims(2) = .FALSE.
	243	CALL MPI_CART_SUB( comm2d, remain_dims, comm1dx, ierr )
	244	CALL MPI_COMM_RANK( comm1dx, myidx, ierr )
	245	!
	246	!-- Transposition from x to y
	247	remain_dims(1) = .FALSE.
	248	remain_dims(2) = .TRUE.
	249	CALL MPI_CART_SUB( comm2d, remain_dims, comm1dy, ierr )
	250	CALL MPI_COMM_RANK( comm1dy, myidy, ierr )
	251
[807]	252	#endif
[1]	253
	254	!
[1003]	255	!-- Calculate array bounds along x-direction for every PE.
[1]	256	ALLOCATE( nxlf(0:pdims(1)-1), nxrf(0:pdims(1)-1), nynf(0:pdims(2)-1), &
[1003]	257	nysf(0:pdims(2)-1) )
[1]	258
[1003]	259	IF ( MOD( nx+1 , pdims(1) ) /= 0 ) THEN
[274]	260	WRITE( message_string, * ) 'x-direction: gridpoint number (',nx+1,') ',&
	261	'is not an& integral divisor of the number ', &
	262	'processors (', pdims(1),')'
[254]	263	CALL message( 'init_pegrid', 'PA0225', 1, 2, 0, 6, 0 )
[1]	264	ELSE
[1003]	265	nnx = ( nx + 1 ) / pdims(1)
[1]	266	IF ( nnx*pdims(1) - ( nx + 1) > nnx ) THEN
[274]	267	WRITE( message_string, * ) 'x-direction: nx does not match the', &
	268	'requirements given by the number of PEs &used', &
	269	'& please use nx = ', nx - ( pdims(1) - ( nnx*pdims(1) &
	270	- ( nx + 1 ) ) ), ' instead of nx =', nx
[254]	271	CALL message( 'init_pegrid', 'PA0226', 1, 2, 0, 6, 0 )
[1]	272	ENDIF
	273	ENDIF
	274
	275	!
	276	!-- Left and right array bounds, number of gridpoints
	277	DO i = 0, pdims(1)-1
	278	nxlf(i) = i * nnx
	279	nxrf(i) = ( i + 1 ) * nnx - 1
	280	ENDDO
	281
	282	!
	283	!-- Calculate array bounds in y-direction for every PE.
[1003]	284	IF ( MOD( ny+1 , pdims(2) ) /= 0 ) THEN
[274]	285	WRITE( message_string, * ) 'y-direction: gridpoint number (',ny+1,') ', &
	286	'is not an& integral divisor of the number of', &
	287	'processors (', pdims(2),')'
[254]	288	CALL message( 'init_pegrid', 'PA0227', 1, 2, 0, 6, 0 )
[1]	289	ELSE
[1003]	290	nny = ( ny + 1 ) / pdims(2)
[1]	291	IF ( nny*pdims(2) - ( ny + 1) > nny ) THEN
[274]	292	WRITE( message_string, * ) 'y-direction: ny does not match the', &
	293	'requirements given by the number of PEs &used ', &
	294	'& please use ny = ', ny - ( pdims(2) - ( nnx*pdims(2) &
[254]	295	- ( ny + 1 ) ) ), ' instead of ny =', ny
	296	CALL message( 'init_pegrid', 'PA0228', 1, 2, 0, 6, 0 )
[1]	297	ENDIF
	298	ENDIF
	299
	300	!
	301	!-- South and north array bounds
	302	DO j = 0, pdims(2)-1
	303	nysf(j) = j * nny
	304	nynf(j) = ( j + 1 ) * nny - 1
	305	ENDDO
	306
	307	!
	308	!-- Local array bounds of the respective PEs
[1003]	309	nxl = nxlf(pcoord(1))
	310	nxr = nxrf(pcoord(1))
	311	nys = nysf(pcoord(2))
	312	nyn = nynf(pcoord(2))
	313	nzb = 0
	314	nzt = nz
	315	nnz = nz
[1]	316
	317	!
[707]	318	!-- Set switches to define if the PE is situated at the border of the virtual
	319	!-- processor grid
	320	IF ( nxl == 0 ) left_border_pe = .TRUE.
	321	IF ( nxr == nx ) right_border_pe = .TRUE.
	322	IF ( nys == 0 ) south_border_pe = .TRUE.
	323	IF ( nyn == ny ) north_border_pe = .TRUE.
	324
	325	!
[1]	326	!-- Calculate array bounds and gridpoint numbers for the transposed arrays
	327	!-- (needed in the pressure solver)
	328	!-- For the transposed arrays, cyclic boundaries as well as top and bottom
	329	!-- boundaries are omitted, because they are obstructive to the transposition
	330
	331	!
	332	!-- 1. transposition z --> x
[1001]	333	!-- This transposition is not neccessary in case of a 1d-decomposition along x
[1304]	334	nys_x = nys
	335	nyn_x = nyn
	336	nny_x = nny
	337	nnz_x = nz / pdims(1)
	338	nzb_x = 1 + myidx * nnz_x
	339	nzt_x = ( myidx + 1 ) * nnz_x
	340	sendrecvcount_zx = nnx * nny * nnz_x
	341
[1001]	342	IF ( pdims(2) /= 1 ) THEN
[1003]	343	IF ( MOD( nz , pdims(1) ) /= 0 ) THEN
[274]	344	WRITE( message_string, * ) 'transposition z --> x:', &
	345	'&nz=',nz,' is not an integral divisior of pdims(1)=', &
	346	pdims(1)
[254]	347	CALL message( 'init_pegrid', 'PA0230', 1, 2, 0, 6, 0 )
[1]	348	ENDIF
	349	ENDIF
	350
	351	!
	352	!-- 2. transposition x --> y
[1003]	353	nnz_y = nnz_x
	354	nzb_y = nzb_x
	355	nzt_y = nzt_x
	356	IF ( MOD( nx+1 , pdims(2) ) /= 0 ) THEN
[274]	357	WRITE( message_string, * ) 'transposition x --> y:', &
	358	'&nx+1=',nx+1,' is not an integral divisor of ',&
	359	'pdims(2)=',pdims(2)
[254]	360	CALL message( 'init_pegrid', 'PA0231', 1, 2, 0, 6, 0 )
[1]	361	ENDIF
[1003]	362	nnx_y = (nx+1) / pdims(2)
[1]	363	nxl_y = myidy * nnx_y
[1003]	364	nxr_y = ( myidy + 1 ) * nnx_y - 1
[1]	365	sendrecvcount_xy = nnx_y * nny_x * nnz_y
	366
	367	!
	368	!-- 3. transposition y --> z (ELSE: x --> y in case of 1D-decomposition
	369	!-- along x)
[1304]	370	nnx_z = nnx_y
	371	nxl_z = nxl_y
	372	nxr_z = nxr_y
	373	nny_z = (ny+1) / pdims(1)
	374	nys_z = myidx * nny_z
	375	nyn_z = ( myidx + 1 ) * nny_z - 1
	376	sendrecvcount_yz = nnx_y * nny_z * nnz_y
	377
[1001]	378	IF ( pdims(2) /= 1 ) THEN
[1]	379	!
	380	!-- y --> z
	381	!-- This transposition is not neccessary in case of a 1d-decomposition
	382	!-- along x, except that the uptream-spline method is switched on
[1003]	383	IF ( MOD( ny+1 , pdims(1) ) /= 0 ) THEN
[274]	384	WRITE( message_string, * ) 'transposition y --> z:', &
	385	'& ny+1=',ny+1,' is not an integral divisor of ',&
	386	'pdims(1)=',pdims(1)
[254]	387	CALL message( 'init_pegrid', 'PA0232', 1, 2, 0, 6, 0 )
[1]	388	ENDIF
	389
	390	ELSE
	391	!
	392	!-- x --> y. This condition must be fulfilled for a 1D-decomposition along x
[1003]	393	IF ( MOD( ny+1 , pdims(1) ) /= 0 ) THEN
[274]	394	WRITE( message_string, * ) 'transposition x --> y:', &
	395	'& ny+1=',ny+1,' is not an integral divisor of ',&
	396	'pdims(1)=',pdims(1)
[254]	397	CALL message( 'init_pegrid', 'PA0233', 1, 2, 0, 6, 0 )
[1]	398	ENDIF
	399
	400	ENDIF
	401
	402	!
	403	!-- Indices for direct transpositions z --> y (used for calculating spectra)
[1353]	404	IF ( dt_dosp /= 9999999.9_wp ) THEN
[1003]	405	IF ( MOD( nz, pdims(2) ) /= 0 ) THEN
[274]	406	WRITE( message_string, * ) 'direct transposition z --> y (needed ', &
	407	'for spectra):& nz=',nz,' is not an integral divisor of ',&
	408	'pdims(2)=',pdims(2)
[254]	409	CALL message( 'init_pegrid', 'PA0234', 1, 2, 0, 6, 0 )
[1]	410	ELSE
[1003]	411	nxl_yd = nxl
	412	nxr_yd = nxr
	413	nzb_yd = 1 + myidy * ( nz / pdims(2) )
	414	nzt_yd = ( myidy + 1 ) * ( nz / pdims(2) )
	415	sendrecvcount_zyd = nnx * nny * ( nz / pdims(2) )
[1]	416	ENDIF
	417	ENDIF
	418
	419	!
	420	!-- Indices for direct transpositions y --> x (they are only possible in case
	421	!-- of a 1d-decomposition along x)
	422	IF ( pdims(2) == 1 ) THEN
[1003]	423	nny_x = nny / pdims(1)
	424	nys_x = myid * nny_x
	425	nyn_x = ( myid + 1 ) * nny_x - 1
	426	nzb_x = 1
	427	nzt_x = nz
	428	sendrecvcount_xy = nnx * nny_x * nz
[1]	429	ENDIF
	430
	431	!
	432	!-- Indices for direct transpositions x --> y (they are only possible in case
	433	!-- of a 1d-decomposition along y)
	434	IF ( pdims(1) == 1 ) THEN
[1003]	435	nnx_y = nnx / pdims(2)
	436	nxl_y = myid * nnx_y
	437	nxr_y = ( myid + 1 ) * nnx_y - 1
	438	nzb_y = 1
	439	nzt_y = nz
	440	sendrecvcount_xy = nnx_y * nny * nz
[1]	441	ENDIF
	442
	443	!
	444	!-- Arrays for storing the array bounds are needed any more
	445	DEALLOCATE( nxlf , nxrf , nynf , nysf )
	446
[807]	447
[809]	448	#if ! defined( __check)
[145]	449	!
	450	!-- Collect index bounds from other PEs (to be written to restart file later)
	451	ALLOCATE( hor_index_bounds(4,0:numprocs-1) )
	452
	453	IF ( myid == 0 ) THEN
	454
	455	hor_index_bounds(1,0) = nxl
	456	hor_index_bounds(2,0) = nxr
	457	hor_index_bounds(3,0) = nys
	458	hor_index_bounds(4,0) = nyn
	459
	460	!
	461	!-- Receive data from all other PEs
	462	DO i = 1, numprocs-1
	463	CALL MPI_RECV( ibuf, 4, MPI_INTEGER, i, MPI_ANY_TAG, comm2d, status, &
	464	ierr )
	465	hor_index_bounds(:,i) = ibuf(1:4)
	466	ENDDO
	467
	468	ELSE
	469	!
	470	!-- Send index bounds to PE0
	471	ibuf(1) = nxl
	472	ibuf(2) = nxr
	473	ibuf(3) = nys
	474	ibuf(4) = nyn
	475	CALL MPI_SEND( ibuf, 4, MPI_INTEGER, 0, myid, comm2d, ierr )
	476
	477	ENDIF
	478
[807]	479	#endif
	480
[1]	481	#if defined( __print )
	482	!
	483	!-- Control output
	484	IF ( myid == 0 ) THEN
	485	PRINT, ' processor topology *'
	486	PRINT*, ' '
	487	PRINT*, 'myid pcoord left right south north idx idy nxl: nxr',&
	488	&' nys: nyn'
	489	PRINT*, '------------------------------------------------------------',&
	490	&'-----------'
	491	WRITE (*,1000) 0, pcoord(1), pcoord(2), pleft, pright, psouth, pnorth, &
	492	myidx, myidy, nxl, nxr, nys, nyn
	493	1000 FORMAT (I4,2X,'(',I3,',',I3,')',3X,I4,2X,I4,3X,I4,2X,I4,2X,I3,1X,I3, &
	494	2(2X,I4,':',I4))
	495
	496	!
[108]	497	!-- Receive data from the other PEs
[1]	498	DO i = 1,numprocs-1
	499	CALL MPI_RECV( ibuf, 12, MPI_INTEGER, i, MPI_ANY_TAG, comm2d, status, &
	500	ierr )
	501	WRITE (*,1000) i, ( ibuf(j) , j = 1,12 )
	502	ENDDO
	503	ELSE
	504
	505	!
	506	!-- Send data to PE0
	507	ibuf(1) = pcoord(1); ibuf(2) = pcoord(2); ibuf(3) = pleft
	508	ibuf(4) = pright; ibuf(5) = psouth; ibuf(6) = pnorth; ibuf(7) = myidx
	509	ibuf(8) = myidy; ibuf(9) = nxl; ibuf(10) = nxr; ibuf(11) = nys
	510	ibuf(12) = nyn
	511	CALL MPI_SEND( ibuf, 12, MPI_INTEGER, 0, myid, comm2d, ierr )
	512	ENDIF
	513	#endif
	514
[809]	515	#if defined( __parallel ) && ! defined( __check)
[102]	516	#if defined( __mpi2 )
	517	!
	518	!-- In case of coupled runs, get the port name on PE0 of the atmosphere model
	519	!-- and pass it to PE0 of the ocean model
	520	IF ( myid == 0 ) THEN
	521
	522	IF ( coupling_mode == 'atmosphere_to_ocean' ) THEN
	523
	524	CALL MPI_OPEN_PORT( MPI_INFO_NULL, port_name, ierr )
[108]	525
[102]	526	CALL MPI_PUBLISH_NAME( 'palm_coupler', MPI_INFO_NULL, port_name, &
	527	ierr )
[108]	528
	529	!
[104]	530	!-- Write a flag file for the ocean model and the other atmosphere
	531	!-- processes.
	532	!-- There seems to be a bug in MPICH2 which causes hanging processes
	533	!-- in case that execution of LOOKUP_NAME is continued too early
	534	!-- (i.e. before the port has been created)
	535	OPEN( 90, FILE='COUPLING_PORT_OPENED', FORM='FORMATTED' )
	536	WRITE ( 90, '(''TRUE'')' )
	537	CLOSE ( 90 )
[102]	538
	539	ELSEIF ( coupling_mode == 'ocean_to_atmosphere' ) THEN
	540
[104]	541	!
	542	!-- Continue only if the atmosphere model has created the port.
	543	!-- There seems to be a bug in MPICH2 which causes hanging processes
	544	!-- in case that execution of LOOKUP_NAME is continued too early
	545	!-- (i.e. before the port has been created)
	546	INQUIRE( FILE='COUPLING_PORT_OPENED', EXIST=found )
	547	DO WHILE ( .NOT. found )
	548	INQUIRE( FILE='COUPLING_PORT_OPENED', EXIST=found )
	549	ENDDO
	550
[102]	551	CALL MPI_LOOKUP_NAME( 'palm_coupler', MPI_INFO_NULL, port_name, ierr )
	552
	553	ENDIF
	554
	555	ENDIF
	556
	557	!
	558	!-- In case of coupled runs, establish the connection between the atmosphere
	559	!-- and the ocean model and define the intercommunicator (comm_inter)
	560	CALL MPI_BARRIER( comm2d, ierr )
	561	IF ( coupling_mode == 'atmosphere_to_ocean' ) THEN
	562
	563	CALL MPI_COMM_ACCEPT( port_name, MPI_INFO_NULL, 0, MPI_COMM_WORLD, &
	564	comm_inter, ierr )
[108]	565	coupling_mode_remote = 'ocean_to_atmosphere'
	566
[102]	567	ELSEIF ( coupling_mode == 'ocean_to_atmosphere' ) THEN
	568
	569	CALL MPI_COMM_CONNECT( port_name, MPI_INFO_NULL, 0, MPI_COMM_WORLD, &
	570	comm_inter, ierr )
[108]	571	coupling_mode_remote = 'atmosphere_to_ocean'
	572
[102]	573	ENDIF
[206]	574	#endif
[102]	575
[667]	576	!
[709]	577	!-- Determine the number of ghost point layers
	578	IF ( scalar_advec == 'ws-scheme' .OR. momentum_advec == 'ws-scheme' ) THEN
[667]	579	nbgp = 3
	580	ELSE
	581	nbgp = 1
[709]	582	ENDIF
[667]	583
[102]	584	!
[709]	585	!-- Create a new MPI derived datatype for the exchange of surface (xy) data,
	586	!-- which is needed for coupled atmosphere-ocean runs.
	587	!-- First, calculate number of grid points of an xy-plane.
[667]	588	ngp_xy = ( nxr - nxl + 1 + 2 * nbgp ) * ( nyn - nys + 1 + 2 * nbgp )
[102]	589	CALL MPI_TYPE_VECTOR( ngp_xy, 1, nzt-nzb+2, MPI_REAL, type_xy, ierr )
	590	CALL MPI_TYPE_COMMIT( type_xy, ierr )
[667]	591
[709]	592	IF ( TRIM( coupling_mode ) /= 'uncoupled' ) THEN
[667]	593
	594	!
	595	!-- Pass the number of grid points of the atmosphere model to
	596	!-- the ocean model and vice versa
	597	IF ( coupling_mode == 'atmosphere_to_ocean' ) THEN
	598
	599	nx_a = nx
	600	ny_a = ny
	601
[709]	602	IF ( myid == 0 ) THEN
	603
	604	CALL MPI_SEND( nx_a, 1, MPI_INTEGER, numprocs, 1, comm_inter, &
	605	ierr )
	606	CALL MPI_SEND( ny_a, 1, MPI_INTEGER, numprocs, 2, comm_inter, &
	607	ierr )
	608	CALL MPI_SEND( pdims, 2, MPI_INTEGER, numprocs, 3, comm_inter, &
	609	ierr )
	610	CALL MPI_RECV( nx_o, 1, MPI_INTEGER, numprocs, 4, comm_inter, &
	611	status, ierr )
	612	CALL MPI_RECV( ny_o, 1, MPI_INTEGER, numprocs, 5, comm_inter, &
	613	status, ierr )
	614	CALL MPI_RECV( pdims_remote, 2, MPI_INTEGER, numprocs, 6, &
[667]	615	comm_inter, status, ierr )
	616	ENDIF
	617
[709]	618	CALL MPI_BCAST( nx_o, 1, MPI_INTEGER, 0, comm2d, ierr )
	619	CALL MPI_BCAST( ny_o, 1, MPI_INTEGER, 0, comm2d, ierr )
	620	CALL MPI_BCAST( pdims_remote, 2, MPI_INTEGER, 0, comm2d, ierr )
[667]	621
	622	ELSEIF ( coupling_mode == 'ocean_to_atmosphere' ) THEN
	623
	624	nx_o = nx
	625	ny_o = ny
	626
	627	IF ( myid == 0 ) THEN
[709]	628
	629	CALL MPI_RECV( nx_a, 1, MPI_INTEGER, 0, 1, comm_inter, status, &
	630	ierr )
	631	CALL MPI_RECV( ny_a, 1, MPI_INTEGER, 0, 2, comm_inter, status, &
	632	ierr )
	633	CALL MPI_RECV( pdims_remote, 2, MPI_INTEGER, 0, 3, comm_inter, &
	634	status, ierr )
	635	CALL MPI_SEND( nx_o, 1, MPI_INTEGER, 0, 4, comm_inter, ierr )
	636	CALL MPI_SEND( ny_o, 1, MPI_INTEGER, 0, 5, comm_inter, ierr )
	637	CALL MPI_SEND( pdims, 2, MPI_INTEGER, 0, 6, comm_inter, ierr )
[667]	638	ENDIF
	639
	640	CALL MPI_BCAST( nx_a, 1, MPI_INTEGER, 0, comm2d, ierr)
	641	CALL MPI_BCAST( ny_a, 1, MPI_INTEGER, 0, comm2d, ierr)
	642	CALL MPI_BCAST( pdims_remote, 2, MPI_INTEGER, 0, comm2d, ierr)
	643
	644	ENDIF
	645
[709]	646	ngp_a = ( nx_a+1 + 2 * nbgp ) * ( ny_a+1 + 2 * nbgp )
	647	ngp_o = ( nx_o+1 + 2 * nbgp ) * ( ny_o+1 + 2 * nbgp )
[667]	648
	649	!
[709]	650	!-- Determine if the horizontal grid and the number of PEs in ocean and
	651	!-- atmosphere is same or not
	652	IF ( nx_o == nx_a .AND. ny_o == ny_a .AND. &
[667]	653	pdims(1) == pdims_remote(1) .AND. pdims(2) == pdims_remote(2) ) &
	654	THEN
	655	coupling_topology = 0
	656	ELSE
	657	coupling_topology = 1
	658	ENDIF
	659
	660	!
	661	!-- Determine the target PEs for the exchange between ocean and
	662	!-- atmosphere (comm2d)
[709]	663	IF ( coupling_topology == 0 ) THEN
	664	!
	665	!-- In case of identical topologies, every atmosphere PE has exactly one
	666	!-- ocean PE counterpart and vice versa
	667	IF ( TRIM( coupling_mode ) == 'atmosphere_to_ocean' ) THEN
[667]	668	target_id = myid + numprocs
	669	ELSE
	670	target_id = myid
	671	ENDIF
	672
	673	ELSE
	674	!
	675	!-- In case of nonequivalent topology in ocean and atmosphere only for
	676	!-- PE0 in ocean and PE0 in atmosphere a target_id is needed, since
[709]	677	!-- data echxchange between ocean and atmosphere will be done only
	678	!-- between these PEs.
	679	IF ( myid == 0 ) THEN
	680
	681	IF ( TRIM( coupling_mode ) == 'atmosphere_to_ocean' ) THEN
[667]	682	target_id = numprocs
	683	ELSE
	684	target_id = 0
	685	ENDIF
[709]	686
[667]	687	ENDIF
[709]	688
[667]	689	ENDIF
	690
	691	ENDIF
	692
	693
[102]	694	#endif
	695
[1]	696	#else
	697
	698	!
	699	!-- Array bounds when running on a single PE (respectively a non-parallel
	700	!-- machine)
[1003]	701	nxl = 0
	702	nxr = nx
	703	nnx = nxr - nxl + 1
	704	nys = 0
	705	nyn = ny
	706	nny = nyn - nys + 1
	707	nzb = 0
	708	nzt = nz
	709	nnz = nz
[1]	710
[145]	711	ALLOCATE( hor_index_bounds(4,0:0) )
	712	hor_index_bounds(1,0) = nxl
	713	hor_index_bounds(2,0) = nxr
	714	hor_index_bounds(3,0) = nys
	715	hor_index_bounds(4,0) = nyn
	716
[1]	717	!
	718	!-- Array bounds for the pressure solver (in the parallel code, these bounds
	719	!-- are the ones for the transposed arrays)
[1003]	720	nys_x = nys
	721	nyn_x = nyn
	722	nzb_x = nzb + 1
	723	nzt_x = nzt
[1]	724
[1003]	725	nxl_y = nxl
	726	nxr_y = nxr
	727	nzb_y = nzb + 1
	728	nzt_y = nzt
[1]	729
[1003]	730	nxl_z = nxl
	731	nxr_z = nxr
	732	nys_z = nys
	733	nyn_z = nyn
[1]	734
	735	#endif
	736
	737	!
	738	!-- Calculate number of grid levels necessary for the multigrid poisson solver
	739	!-- as well as the gridpoint indices on each level
	740	IF ( psolver == 'multigrid' ) THEN
	741
	742	!
	743	!-- First calculate number of possible grid levels for the subdomains
	744	mg_levels_x = 1
	745	mg_levels_y = 1
	746	mg_levels_z = 1
	747
	748	i = nnx
	749	DO WHILE ( MOD( i, 2 ) == 0 .AND. i /= 2 )
	750	i = i / 2
	751	mg_levels_x = mg_levels_x + 1
	752	ENDDO
	753
	754	j = nny
	755	DO WHILE ( MOD( j, 2 ) == 0 .AND. j /= 2 )
	756	j = j / 2
	757	mg_levels_y = mg_levels_y + 1
	758	ENDDO
	759
[181]	760	k = nz ! do not use nnz because it might be > nz due to transposition
	761	! requirements
[1]	762	DO WHILE ( MOD( k, 2 ) == 0 .AND. k /= 2 )
	763	k = k / 2
	764	mg_levels_z = mg_levels_z + 1
	765	ENDDO
	766
	767	maximum_grid_level = MIN( mg_levels_x, mg_levels_y, mg_levels_z )
	768
	769	!
	770	!-- Find out, if the total domain allows more levels. These additional
[709]	771	!-- levels are identically processed on all PEs.
[197]	772	IF ( numprocs > 1 .AND. mg_switch_to_pe0_level /= -1 ) THEN
[709]	773
[1]	774	IF ( mg_levels_z > MIN( mg_levels_x, mg_levels_y ) ) THEN
[709]	775
[1]	776	mg_switch_to_pe0_level_l = maximum_grid_level
	777
	778	mg_levels_x = 1
	779	mg_levels_y = 1
	780
	781	i = nx+1
	782	DO WHILE ( MOD( i, 2 ) == 0 .AND. i /= 2 )
	783	i = i / 2
	784	mg_levels_x = mg_levels_x + 1
	785	ENDDO
	786
	787	j = ny+1
	788	DO WHILE ( MOD( j, 2 ) == 0 .AND. j /= 2 )
	789	j = j / 2
	790	mg_levels_y = mg_levels_y + 1
	791	ENDDO
	792
	793	maximum_grid_level_l = MIN( mg_levels_x, mg_levels_y, mg_levels_z )
	794
	795	IF ( maximum_grid_level_l > mg_switch_to_pe0_level_l ) THEN
	796	mg_switch_to_pe0_level_l = maximum_grid_level_l - &
	797	mg_switch_to_pe0_level_l + 1
	798	ELSE
	799	mg_switch_to_pe0_level_l = 0
	800	ENDIF
[709]	801
[1]	802	ELSE
	803	mg_switch_to_pe0_level_l = 0
	804	maximum_grid_level_l = maximum_grid_level
[709]	805
[1]	806	ENDIF
	807
	808	!
	809	!-- Use switch level calculated above only if it is not pre-defined
	810	!-- by user
	811	IF ( mg_switch_to_pe0_level == 0 ) THEN
	812	IF ( mg_switch_to_pe0_level_l /= 0 ) THEN
	813	mg_switch_to_pe0_level = mg_switch_to_pe0_level_l
	814	maximum_grid_level = maximum_grid_level_l
	815	ENDIF
	816
	817	ELSE
	818	!
	819	!-- Check pre-defined value and reset to default, if neccessary
	820	IF ( mg_switch_to_pe0_level < mg_switch_to_pe0_level_l .OR. &
	821	mg_switch_to_pe0_level >= maximum_grid_level_l ) THEN
[254]	822	message_string = 'mg_switch_to_pe0_level ' // &
	823	'out of range and reset to default (=0)'
	824	CALL message( 'init_pegrid', 'PA0235', 0, 1, 0, 6, 0 )
[1]	825	mg_switch_to_pe0_level = 0
	826	ELSE
	827	!
	828	!-- Use the largest number of possible levels anyway and recalculate
	829	!-- the switch level to this largest number of possible values
	830	maximum_grid_level = maximum_grid_level_l
	831
	832	ENDIF
[709]	833
[1]	834	ENDIF
	835
	836	ENDIF
	837
[1056]	838	ALLOCATE( grid_level_count(maximum_grid_level), &
	839	nxl_mg(0:maximum_grid_level), nxr_mg(0:maximum_grid_level), &
	840	nyn_mg(0:maximum_grid_level), nys_mg(0:maximum_grid_level), &
	841	nzt_mg(0:maximum_grid_level) )
[1]	842
	843	grid_level_count = 0
[1056]	844	!
	845	!-- Index zero required as dummy due to definition of arrays f2 and p2 in
	846	!-- recursive subroutine next_mg_level
	847	nxl_mg(0) = 0; nxr_mg(0) = 0; nyn_mg(0) = 0; nys_mg(0) = 0; nzt_mg(0) = 0
[778]	848
[1]	849	nxl_l = nxl; nxr_l = nxr; nys_l = nys; nyn_l = nyn; nzt_l = nzt
	850
	851	DO i = maximum_grid_level, 1 , -1
	852
	853	IF ( i == mg_switch_to_pe0_level ) THEN
[809]	854	#if defined( __parallel ) && ! defined( __check )
[1]	855	!
	856	!-- Save the grid size of the subdomain at the switch level, because
	857	!-- it is needed in poismg.
	858	ind(1) = nxl_l; ind(2) = nxr_l
	859	ind(3) = nys_l; ind(4) = nyn_l
	860	ind(5) = nzt_l
	861	ALLOCATE( ind_all(5*numprocs), mg_loc_ind(5,0:numprocs-1) )
	862	CALL MPI_ALLGATHER( ind, 5, MPI_INTEGER, ind_all, 5, &
	863	MPI_INTEGER, comm2d, ierr )
	864	DO j = 0, numprocs-1
	865	DO k = 1, 5
	866	mg_loc_ind(k,j) = ind_all(k+j*5)
	867	ENDDO
	868	ENDDO
	869	DEALLOCATE( ind_all )
	870	!
[709]	871	!-- Calculate the grid size of the total domain
[1]	872	nxr_l = ( nxr_l-nxl_l+1 ) * pdims(1) - 1
	873	nxl_l = 0
	874	nyn_l = ( nyn_l-nys_l+1 ) * pdims(2) - 1
	875	nys_l = 0
	876	!
	877	!-- The size of this gathered array must not be larger than the
	878	!-- array tend, which is used in the multigrid scheme as a temporary
[778]	879	!-- array. Therefore the subdomain size of an PE is calculated and
	880	!-- the size of the gathered grid. These values are used in
	881	!-- routines pres and poismg
	882	subdomain_size = ( nxr - nxl + 2 * nbgp + 1 ) * &
	883	( nyn - nys + 2 * nbgp + 1 ) * ( nzt - nzb + 2 )
[1]	884	gathered_size = ( nxr_l - nxl_l + 3 ) * ( nyn_l - nys_l + 3 ) * &
	885	( nzt_l - nzb + 2 )
	886
[809]	887	#elif ! defined ( __parallel )
[254]	888	message_string = 'multigrid gather/scatter impossible ' // &
[1]	889	'in non parallel mode'
[254]	890	CALL message( 'init_pegrid', 'PA0237', 1, 2, 0, 6, 0 )
[1]	891	#endif
	892	ENDIF
	893
	894	nxl_mg(i) = nxl_l
	895	nxr_mg(i) = nxr_l
	896	nys_mg(i) = nys_l
	897	nyn_mg(i) = nyn_l
	898	nzt_mg(i) = nzt_l
	899
	900	nxl_l = nxl_l / 2
	901	nxr_l = nxr_l / 2
	902	nys_l = nys_l / 2
	903	nyn_l = nyn_l / 2
	904	nzt_l = nzt_l / 2
[778]	905
[1]	906	ENDDO
	907
[780]	908	!
	909	!-- Temporary problem: Currently calculation of maxerror iin routine poismg crashes
	910	!-- if grid data are collected on PE0 already on the finest grid level.
	911	!-- To be solved later.
	912	IF ( maximum_grid_level == mg_switch_to_pe0_level ) THEN
	913	message_string = 'grid coarsening on subdomain level cannot be performed'
	914	CALL message( 'poismg', 'PA0236', 1, 2, 0, 6, 0 )
	915	ENDIF
	916
[1]	917	ELSE
	918
[667]	919	maximum_grid_level = 0
[1]	920
	921	ENDIF
	922
[722]	923	!
	924	!-- Default level 0 tells exchange_horiz that all ghost planes have to be
	925	!-- exchanged. grid_level is adjusted in poismg, where only one ghost plane
	926	!-- is required.
	927	grid_level = 0
[1]	928
[809]	929	#if defined( __parallel ) && ! defined ( __check )
[1]	930	!
	931	!-- Gridpoint number for the exchange of ghost points (y-line for 2D-arrays)
[667]	932	ngp_y = nyn - nys + 1 + 2 * nbgp
[1]	933
	934	!
[709]	935	!-- Define new MPI derived datatypes for the exchange of ghost points in
	936	!-- x- and y-direction for 2D-arrays (line)
	937	CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp, ngp_y, MPI_REAL, type_x, &
	938	ierr )
[1]	939	CALL MPI_TYPE_COMMIT( type_x, ierr )
[709]	940	CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp, ngp_y, MPI_INTEGER, &
	941	type_x_int, ierr )
[1]	942	CALL MPI_TYPE_COMMIT( type_x_int, ierr )
	943
[667]	944	CALL MPI_TYPE_VECTOR( nbgp, ngp_y, ngp_y, MPI_REAL, type_y, ierr )
	945	CALL MPI_TYPE_COMMIT( type_y, ierr )
	946	CALL MPI_TYPE_VECTOR( nbgp, ngp_y, ngp_y, MPI_INTEGER, type_y_int, ierr )
	947	CALL MPI_TYPE_COMMIT( type_y_int, ierr )
	948
	949
[1]	950	!
	951	!-- Calculate gridpoint numbers for the exchange of ghost points along x
	952	!-- (yz-plane for 3D-arrays) and define MPI derived data type(s) for the
	953	!-- exchange of ghost points in y-direction (xz-plane).
	954	!-- Do these calculations for the model grid and (if necessary) also
	955	!-- for the coarser grid levels used in the multigrid method
[667]	956	ALLOCATE ( ngp_yz(0:maximum_grid_level), type_xz(0:maximum_grid_level),&
	957	type_yz(0:maximum_grid_level) )
[1]	958
	959	nxl_l = nxl; nxr_l = nxr; nys_l = nys; nyn_l = nyn; nzb_l = nzb; nzt_l = nzt
[709]	960
[667]	961	!
	962	!-- Discern between the model grid, which needs nbgp ghost points and
	963	!-- grid levels for the multigrid scheme. In the latter case only one
	964	!-- ghost point is necessary.
[709]	965	!-- First definition of MPI-datatypes for exchange of ghost layers on normal
[667]	966	!-- grid. The following loop is needed for data exchange in poismg.f90.
	967	!
	968	!-- Determine number of grid points of yz-layer for exchange
	969	ngp_yz(0) = (nzt - nzb + 2) * (nyn - nys + 1 + 2 * nbgp)
[709]	970
[667]	971	!
[709]	972	!-- Define an MPI-datatype for the exchange of left/right boundaries.
	973	!-- Although data are contiguous in physical memory (which does not
	974	!-- necessarily require an MPI-derived datatype), the data exchange between
	975	!-- left and right PE's using the MPI-derived type is 10% faster than without.
[667]	976	CALL MPI_TYPE_VECTOR( nxr-nxl+1+2nbgp, nbgp(nzt-nzb+2), ngp_yz(0), &
[709]	977	MPI_REAL, type_xz(0), ierr )
[667]	978	CALL MPI_TYPE_COMMIT( type_xz(0), ierr )
[1]	979
[709]	980	CALL MPI_TYPE_VECTOR( nbgp, ngp_yz(0), ngp_yz(0), MPI_REAL, type_yz(0), &
	981	ierr )
[667]	982	CALL MPI_TYPE_COMMIT( type_yz(0), ierr )
[709]	983
[667]	984	!
[709]	985	!-- Definition of MPI-datatypes for multigrid method (coarser level grids)
[667]	986	IF ( psolver == 'multigrid' ) THEN
	987	!
[709]	988	!-- Definition of MPI-datatyoe as above, but only 1 ghost level is used
	989	DO i = maximum_grid_level, 1 , -1
	990
[667]	991	ngp_yz(i) = (nzt_l - nzb_l + 2) * (nyn_l - nys_l + 3)
	992
	993	CALL MPI_TYPE_VECTOR( nxr_l-nxl_l+3, nzt_l-nzb_l+2, ngp_yz(i), &
[709]	994	MPI_REAL, type_xz(i), ierr )
[667]	995	CALL MPI_TYPE_COMMIT( type_xz(i), ierr )
[1]	996
[709]	997	CALL MPI_TYPE_VECTOR( 1, ngp_yz(i), ngp_yz(i), MPI_REAL, type_yz(i), &
	998	ierr )
[667]	999	CALL MPI_TYPE_COMMIT( type_yz(i), ierr )
	1000
	1001	nxl_l = nxl_l / 2
	1002	nxr_l = nxr_l / 2
	1003	nys_l = nys_l / 2
	1004	nyn_l = nyn_l / 2
	1005	nzt_l = nzt_l / 2
[709]	1006
[667]	1007	ENDDO
[709]	1008
	1009	ENDIF
[1]	1010	#endif
	1011
[809]	1012	#if defined( __parallel ) && ! defined ( __check )
[1]	1013	!
	1014	!-- Setting of flags for inflow/outflow conditions in case of non-cyclic
[106]	1015	!-- horizontal boundary conditions.
[1]	1016	IF ( pleft == MPI_PROC_NULL ) THEN
[1159]	1017	IF ( bc_lr == 'dirichlet/radiation' ) THEN
[1]	1018	inflow_l = .TRUE.
[1159]	1019	ELSEIF ( bc_lr == 'radiation/dirichlet' ) THEN
[1]	1020	outflow_l = .TRUE.
	1021	ENDIF
	1022	ENDIF
	1023
	1024	IF ( pright == MPI_PROC_NULL ) THEN
[1159]	1025	IF ( bc_lr == 'dirichlet/radiation' ) THEN
[1]	1026	outflow_r = .TRUE.
[1159]	1027	ELSEIF ( bc_lr == 'radiation/dirichlet' ) THEN
[1]	1028	inflow_r = .TRUE.
	1029	ENDIF
	1030	ENDIF
	1031
	1032	IF ( psouth == MPI_PROC_NULL ) THEN
[1159]	1033	IF ( bc_ns == 'dirichlet/radiation' ) THEN
[1]	1034	outflow_s = .TRUE.
[1159]	1035	ELSEIF ( bc_ns == 'radiation/dirichlet' ) THEN
[1]	1036	inflow_s = .TRUE.
	1037	ENDIF
	1038	ENDIF
	1039
	1040	IF ( pnorth == MPI_PROC_NULL ) THEN
[1159]	1041	IF ( bc_ns == 'dirichlet/radiation' ) THEN
[1]	1042	inflow_n = .TRUE.
[1159]	1043	ELSEIF ( bc_ns == 'radiation/dirichlet' ) THEN
[1]	1044	outflow_n = .TRUE.
	1045	ENDIF
	1046	ENDIF
	1047
[151]	1048	!
	1049	!-- Broadcast the id of the inflow PE
	1050	IF ( inflow_l ) THEN
[163]	1051	id_inflow_l = myidx
[151]	1052	ELSE
	1053	id_inflow_l = 0
	1054	ENDIF
[622]	1055	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[151]	1056	CALL MPI_ALLREDUCE( id_inflow_l, id_inflow, 1, MPI_INTEGER, MPI_SUM, &
	1057	comm1dx, ierr )
	1058
[163]	1059	!
	1060	!-- Broadcast the id of the recycling plane
	1061	!-- WARNING: needs to be adjusted in case of inflows other than from left side!
[1139]	1062	IF ( NINT( recycling_width / dx ) >= nxl .AND. &
	1063	NINT( recycling_width / dx ) <= nxr ) THEN
[163]	1064	id_recycling_l = myidx
	1065	ELSE
	1066	id_recycling_l = 0
	1067	ENDIF
[622]	1068	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[163]	1069	CALL MPI_ALLREDUCE( id_recycling_l, id_recycling, 1, MPI_INTEGER, MPI_SUM, &
	1070	comm1dx, ierr )
	1071
[809]	1072	#elif ! defined ( __parallel )
[1159]	1073	IF ( bc_lr == 'dirichlet/radiation' ) THEN
[1]	1074	inflow_l = .TRUE.
	1075	outflow_r = .TRUE.
[1159]	1076	ELSEIF ( bc_lr == 'radiation/dirichlet' ) THEN
[1]	1077	outflow_l = .TRUE.
	1078	inflow_r = .TRUE.
	1079	ENDIF
	1080
[1159]	1081	IF ( bc_ns == 'dirichlet/radiation' ) THEN
[1]	1082	inflow_n = .TRUE.
	1083	outflow_s = .TRUE.
[1159]	1084	ELSEIF ( bc_ns == 'radiation/dirichlet' ) THEN
[1]	1085	outflow_n = .TRUE.
	1086	inflow_s = .TRUE.
	1087	ENDIF
	1088	#endif
[807]	1089
[106]	1090	!
[978]	1091	!-- At the inflow or outflow, u or v, respectively, have to be calculated for
	1092	!-- one more grid point.
	1093	IF ( inflow_l .OR. outflow_l ) THEN
[106]	1094	nxlu = nxl + 1
	1095	ELSE
	1096	nxlu = nxl
	1097	ENDIF
[978]	1098	IF ( inflow_s .OR. outflow_s ) THEN
[106]	1099	nysv = nys + 1
	1100	ELSE
	1101	nysv = nys
	1102	ENDIF
[1]	1103
[114]	1104	!
	1105	!-- Allocate wall flag arrays used in the multigrid solver
	1106	IF ( psolver == 'multigrid' ) THEN
	1107
	1108	DO i = maximum_grid_level, 1, -1
	1109
	1110	SELECT CASE ( i )
	1111
	1112	CASE ( 1 )
	1113	ALLOCATE( wall_flags_1(nzb:nzt_mg(i)+1, &
	1114	nys_mg(i)-1:nyn_mg(i)+1, &
	1115	nxl_mg(i)-1:nxr_mg(i)+1) )
	1116
	1117	CASE ( 2 )
	1118	ALLOCATE( wall_flags_2(nzb:nzt_mg(i)+1, &
	1119	nys_mg(i)-1:nyn_mg(i)+1, &
	1120	nxl_mg(i)-1:nxr_mg(i)+1) )
	1121
	1122	CASE ( 3 )
	1123	ALLOCATE( wall_flags_3(nzb:nzt_mg(i)+1, &
	1124	nys_mg(i)-1:nyn_mg(i)+1, &
	1125	nxl_mg(i)-1:nxr_mg(i)+1) )
	1126
	1127	CASE ( 4 )
	1128	ALLOCATE( wall_flags_4(nzb:nzt_mg(i)+1, &
	1129	nys_mg(i)-1:nyn_mg(i)+1, &
	1130	nxl_mg(i)-1:nxr_mg(i)+1) )
	1131
	1132	CASE ( 5 )
	1133	ALLOCATE( wall_flags_5(nzb:nzt_mg(i)+1, &
	1134	nys_mg(i)-1:nyn_mg(i)+1, &
	1135	nxl_mg(i)-1:nxr_mg(i)+1) )
	1136
	1137	CASE ( 6 )
	1138	ALLOCATE( wall_flags_6(nzb:nzt_mg(i)+1, &
	1139	nys_mg(i)-1:nyn_mg(i)+1, &
	1140	nxl_mg(i)-1:nxr_mg(i)+1) )
	1141
	1142	CASE ( 7 )
	1143	ALLOCATE( wall_flags_7(nzb:nzt_mg(i)+1, &
	1144	nys_mg(i)-1:nyn_mg(i)+1, &
	1145	nxl_mg(i)-1:nxr_mg(i)+1) )
	1146
	1147	CASE ( 8 )
	1148	ALLOCATE( wall_flags_8(nzb:nzt_mg(i)+1, &
	1149	nys_mg(i)-1:nyn_mg(i)+1, &
	1150	nxl_mg(i)-1:nxr_mg(i)+1) )
	1151
	1152	CASE ( 9 )
	1153	ALLOCATE( wall_flags_9(nzb:nzt_mg(i)+1, &
	1154	nys_mg(i)-1:nyn_mg(i)+1, &
	1155	nxl_mg(i)-1:nxr_mg(i)+1) )
	1156
	1157	CASE ( 10 )
	1158	ALLOCATE( wall_flags_10(nzb:nzt_mg(i)+1, &
	1159	nys_mg(i)-1:nyn_mg(i)+1, &
	1160	nxl_mg(i)-1:nxr_mg(i)+1) )
	1161
	1162	CASE DEFAULT
[254]	1163	message_string = 'more than 10 multigrid levels'
	1164	CALL message( 'init_pegrid', 'PA0238', 1, 2, 0, 6, 0 )
[114]	1165
	1166	END SELECT
	1167
	1168	ENDDO
	1169
	1170	ENDIF
	1171
[759]	1172	!
	1173	!-- Calculate the number of groups into which parallel I/O is split.
	1174	!-- The default for files which are opened by all PEs (or where each
	1175	!-- PE opens his own independent file) is, that all PEs are doing input/output
	1176	!-- in parallel at the same time. This might cause performance or even more
	1177	!-- severe problems depending on the configuration of the underlying file
	1178	!-- system.
	1179	!-- First, set the default:
	1180	IF ( maximum_parallel_io_streams == -1 .OR. &
	1181	maximum_parallel_io_streams > numprocs ) THEN
	1182	maximum_parallel_io_streams = numprocs
	1183	ENDIF
	1184
	1185	!
	1186	!-- Now calculate the number of io_blocks and the io_group to which the
	1187	!-- respective PE belongs. I/O of the groups is done in serial, but in parallel
	1188	!-- for all PEs belonging to the same group. A preliminary setting with myid
	1189	!-- based on MPI_COMM_WORLD has been done in parin.
	1190	io_blocks = numprocs / maximum_parallel_io_streams
	1191	io_group = MOD( myid+1, io_blocks )
	1192
	1193
[1]	1194	END SUBROUTINE init_pegrid

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |