Home

Context Navigation

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

Last change on this file since 4286 was 4264, checked in by scharf, 5 years ago
corrected error message string for incompatible combinations of nx, ny or nz with pdims(1) or pdims(2)
Property svn:keywords set to `Id`
File size: 50.6 KB

Rev	Line
[1682]	1	!> @file init_pegrid.f90
[2000]	2	!------------------------------------------------------------------------------!
[2696]	3	! This file is part of the PALM model system.
[1036]	4	!
[2000]	5	! PALM is free software: you can redistribute it and/or modify it under the
	6	! terms of the GNU General Public License as published by the Free Software
	7	! Foundation, either version 3 of the License, or (at your option) any later
	8	! version.
[1036]	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	!
[3655]	17	! Copyright 1997-2019 Leibniz Universitaet Hannover
[2000]	18	!------------------------------------------------------------------------------!
[1036]	19	!
[254]	20	! Current revisions:
[1322]	21	! ------------------
[2201]	22	!
[3553]	23	!
[1321]	24	! Former revisions:
	25	! -----------------
	26	! $Id: init_pegrid.f90 4264 2019-10-15 16:00:23Z resler $
[4264]	27	! corrected error message string
	28	!
	29	! 4241 2019-09-27 06:32:47Z raasch
[4241]	30	! Check added to ensure that subdomain grid has at least the size as given by the number
	31	! of ghost points
	32	!
	33	! 4182 2019-08-22 15:20:23Z scharf
[4182]	34	! Corrected "Former revisions" section
	35	!
	36	! 4045 2019-06-21 10:58:47Z raasch
[4045]	37	! bugfix: kind attribute added to nint function to allow for large integers which may appear in
	38	! case of default recycling width and small grid spacings
	39	!
	40	! 3999 2019-05-23 16:09:37Z suehring
[3999]	41	! Spend 3 ghost points also in case of pw-scheme when nesting is applied
	42	!
	43	! 3897 2019-04-15 11:51:14Z suehring
[3897]	44	! Minor revision of multigrid check; give warning instead of an abort.
	45	!
	46	! 3890 2019-04-12 15:59:20Z suehring
[3890]	47	! Check if grid coarsening is possible on subdomain, in order to avoid that
	48	! multigrid approach effectively reduces to a Gauss-Seidel scheme.
	49	!
	50	! 3885 2019-04-11 11:29:34Z kanani
[3885]	51	! Changes related to global restructuring of location messages and introduction
	52	! of additional debug messages
	53	!
	54	! 3884 2019-04-10 13:31:55Z Giersch
[3884]	55	! id_recycling is only calculated in case of tubulent inflow
	56	!
	57	! 3761 2019-02-25 15:31:42Z raasch
[3761]	58	! unused variable removed
	59	!
	60	! 3655 2019-01-07 16:51:22Z knoop
[3553]	61	! variables documented
[2259]	62	!
[4182]	63	! Revision 1.1 1997/07/24 11:15:09 raasch
	64	! Initial revision
	65	!
	66	!
[1]	67	! Description:
	68	! ------------
[1682]	69	!> Determination of the virtual processor topology (if not prescribed by the
	70	!> user)and computation of the grid point number and array bounds of the local
	71	!> domains.
[2696]	72	!> @todo: remove MPI-data types for 2D exchange on coarse multigrid level (not
	73	!> used any more)
[1]	74	!------------------------------------------------------------------------------!
[1682]	75	SUBROUTINE init_pegrid
	76
[1]	77
[1320]	78	USE control_parameters, &
[3182]	79	ONLY: bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, bc_dirichlet_s, &
	80	bc_lr, bc_ns, bc_radiation_l, bc_radiation_n, bc_radiation_r, &
[3241]	81	bc_radiation_s, coupling_mode, coupling_topology, gathered_size,&
	82	grid_level, grid_level_count, maximum_grid_level, &
[3761]	83	message_string, mg_switch_to_pe0_level, momentum_advec, &
[3182]	84	psolver, outflow_source_plane, recycling_width, scalar_advec, &
[3884]	85	subdomain_size, turbulent_inflow, turbulent_outflow, y_shift
[1]	86
[1320]	87	USE grid_variables, &
	88	ONLY: dx
	89
	90	USE indices, &
	91	ONLY: mg_loc_ind, nbgp, nnx, nny, nnz, nx, nx_a, nx_o, nxl, nxl_mg, &
	92	nxlu, nxr, nxr_mg, ny, ny_a, ny_o, nyn, nyn_mg, nys, nys_mg, &
	93	nysv, nz, nzb, nzt, nzt_mg, wall_flags_1, wall_flags_2, &
	94	wall_flags_3, wall_flags_4, wall_flags_5, wall_flags_6, &
	95	wall_flags_7, wall_flags_8, wall_flags_9, wall_flags_10
[1]	96
[1320]	97	USE kinds
	98
	99	USE pegrid
[3999]	100
	101	USE pmc_interface, &
	102	ONLY: nested_run
[2238]	103
[1833]	104	USE spectra_mod, &
[1922]	105	ONLY: calculate_spectra, dt_dosp
[1833]	106
[2259]	107	USE synthetic_turbulence_generator_mod, &
[2938]	108	ONLY: id_stg_left, id_stg_north, id_stg_right, id_stg_south, &
	109	use_syn_turb_gen
[2259]	110
[1320]	111	USE transpose_indices, &
	112	ONLY: nxl_y, nxl_yd, nxl_z, nxr_y, nxr_yd, nxr_z, nyn_x, nyn_z, nys_x,&
	113	nys_z, nzb_x, nzb_y, nzb_yd, nzt_x, nzt_yd, nzt_y
[667]	114
[2365]	115	USE vertical_nesting_mod, &
	116	ONLY: vnested, vnest_init_pegrid_domain, vnest_init_pegrid_rank
	117
[1]	118	IMPLICIT NONE
	119
[3552]	120	INTEGER(iwp) :: i !< running index over number of processors or number of multigrid level
	121	INTEGER(iwp) :: id_inflow_l !< ID indicating processors located at the left inflow boundary
[2050]	122	INTEGER(iwp) :: id_outflow_l !< local value of id_outflow
	123	INTEGER(iwp) :: id_outflow_source_l !< local value of id_outflow_source
[3552]	124	INTEGER(iwp) :: id_recycling_l !< ID indicating processors located at the recycling plane
[2938]	125	INTEGER(iwp) :: id_stg_left_l !< left lateral boundary local core id in case of turbulence generator
	126	INTEGER(iwp) :: id_stg_north_l !< north lateral boundary local core id in case of turbulence generator
	127	INTEGER(iwp) :: id_stg_right_l !< right lateral boundary local core id in case of turbulence generator
	128	INTEGER(iwp) :: id_stg_south_l !< south lateral boundary local core id in case of turbulence generator
[3552]	129	INTEGER(iwp) :: ind(5) !< array containing the subdomain bounds
	130	INTEGER(iwp) :: j !< running index, used for various loops
	131	INTEGER(iwp) :: k !< number of vertical grid points in different multigrid level
	132	INTEGER(iwp) :: maximum_grid_level_l !< maximum number of grid level without switching to PE 0
	133	INTEGER(iwp) :: mg_levels_x !< maximum number of grid level allowed along x-direction
	134	INTEGER(iwp) :: mg_levels_y !< maximum number of grid level allowed along y-direction
	135	INTEGER(iwp) :: mg_levels_z !< maximum number of grid level allowed along z-direction
	136	INTEGER(iwp) :: mg_switch_to_pe0_level_l !< maximum number of grid level with switching to PE 0
	137	INTEGER(iwp) :: nnx_y !< quotient of number of grid points along x-direction and number of PEs used along y-direction
	138	INTEGER(iwp) :: nny_x !< quotient of number of grid points along y-direction and number of PEs used along x-direction
	139	INTEGER(iwp) :: nny_z !< quotient of number of grid points along y-direction and number of PEs used along x-direction
	140	INTEGER(iwp) :: nnz_x !< quotient of number of grid points along z-direction and number of PEs used along x-direction
	141	INTEGER(iwp) :: nnz_y !< quotient of number of grid points along z-direction and number of PEs used along x-direction
	142	INTEGER(iwp) :: numproc_sqr !< square root of the number of processors
	143	INTEGER(iwp) :: nxl_l !< lower index bound along x-direction on subdomain and different multigrid level
	144	INTEGER(iwp) :: nxr_l !< upper index bound along x-direction on subdomain and different multigrid level
	145	INTEGER(iwp) :: nyn_l !< lower index bound along y-direction on subdomain and different multigrid level
	146	INTEGER(iwp) :: nys_l !< upper index bound along y-direction on subdomain and different multigrid level
	147	INTEGER(iwp) :: nzb_l !< lower index bound along z-direction on subdomain and different multigrid level
	148	INTEGER(iwp) :: nzt_l !< upper index bound along z-direction on subdomain and different multigrid level
	149	!$ INTEGER(iwp) :: omp_get_num_threads !< number of OpenMP threads
[1]	150
[3552]	151	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: ind_all !< dummy array containing index bounds on subdomain, used for gathering
	152	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nxlf !< lower index bound allong x-direction for every PE
	153	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nxrf !< upper index bound allong x-direction for every PE
	154	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nynf !< lower index bound allong y-direction for every PE
	155	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nysf !< lower index bound allong y-direction for every PE
[1]	156
[3552]	157	INTEGER(iwp), DIMENSION(2) :: pdims_remote !< number of PEs used for coupled model (only in atmospher-ocean coupling)
[2372]	158	INTEGER(iwp) :: lcoord(2) !< PE coordinates of left neighbor along x and y
	159	INTEGER(iwp) :: rcoord(2) !< PE coordinates of right neighbor along x and y
[667]	160
[1]	161	!
	162	!-- Get the number of OpenMP threads
	163	!$OMP PARALLEL
	164	!$ threads_per_task = omp_get_num_threads()
	165	!$OMP END PARALLEL
	166
	167
	168	#if defined( __parallel )
[667]	169
[3885]	170	CALL location_message( 'creating virtual PE grids + MPI derived data types', 'start' )
[1764]	171
[1]	172	!
[2177]	173	!-- Determine the processor topology or check it, if prescribed by the user
[1779]	174	IF ( npex == -1 .AND. npey == -1 ) THEN
[1]	175
	176	!
[2177]	177	!-- Automatic determination of the topology
[1779]	178	numproc_sqr = SQRT( REAL( numprocs, KIND=wp ) )
	179	pdims(1) = MAX( numproc_sqr , 1 )
[2180]	180	DO WHILE ( MOD( numprocs , pdims(1) ) /= 0 )
[1779]	181	pdims(1) = pdims(1) - 1
	182	ENDDO
[2180]	183	pdims(2) = numprocs / pdims(1)
[1]	184
[1779]	185	ELSEIF ( npex /= -1 .AND. npey /= -1 ) THEN
[1]	186
	187	!
[1779]	188	!-- Prescribed by user. Number of processors on the prescribed topology
	189	!-- must be equal to the number of PEs available to the job
	190	IF ( ( npex * npey ) /= numprocs ) THEN
[3045]	191	WRITE( message_string, * ) 'number of PEs of the prescribed ', &
[3046]	192	'topology (', npex*npey,') does not match & the number of ', &
[1779]	193	'PEs available to the job (', numprocs, ')'
	194	CALL message( 'init_pegrid', 'PA0221', 1, 2, 0, 6, 0 )
	195	ENDIF
	196	pdims(1) = npex
	197	pdims(2) = npey
[1]	198
[1779]	199	ELSE
[1]	200	!
[1779]	201	!-- If the processor topology is prescribed by the user, the number of
	202	!-- PEs must be given in both directions
[3045]	203	message_string = 'if the processor topology is prescribed by th' // &
[3046]	204	'e user & both values of "npex" and "npey" must be given' // &
[1779]	205	' in the &NAMELIST-parameter file'
	206	CALL message( 'init_pegrid', 'PA0222', 1, 2, 0, 6, 0 )
[1]	207
	208	ENDIF
	209
	210	!
	211	!-- If necessary, set horizontal boundary conditions to non-cyclic
[722]	212	IF ( bc_lr /= 'cyclic' ) cyclic(1) = .FALSE.
	213	IF ( bc_ns /= 'cyclic' ) cyclic(2) = .FALSE.
[1]	214
[807]	215
[1]	216	!
	217	!-- Create the virtual processor grid
	218	CALL MPI_CART_CREATE( comm_palm, ndim, pdims, cyclic, reorder, &
	219	comm2d, ierr )
	220	CALL MPI_COMM_RANK( comm2d, myid, ierr )
[1468]	221	WRITE (myid_char,'(''_'',I6.6)') myid
[1]	222
	223	CALL MPI_CART_COORDS( comm2d, myid, ndim, pcoord, ierr )
	224	CALL MPI_CART_SHIFT( comm2d, 0, 1, pleft, pright, ierr )
	225	CALL MPI_CART_SHIFT( comm2d, 1, 1, psouth, pnorth, ierr )
	226	!
[2372]	227	!-- In case of cyclic boundary conditions, a y-shift at the boundaries in
	228	!-- x-direction can be introduced via parameter y_shift. The shift is done
	229	!-- by modifying the processor grid in such a way that processors located
	230	!-- at the x-boundary communicate across it to processors with y-coordinate
	231	!-- shifted by y_shift relative to their own. This feature can not be used
	232	!-- in combination with an fft pressure solver. It has been implemented to
	233	!-- counter the effect of streak structures in case of cyclic boundary
	234	!-- conditions. For a description of these see Munters
	235	!-- (2016; dx.doi.org/10.1063/1.4941912)
	236	!--
	237	!-- Get coordinates of left and right neighbor on PE grid
	238	IF ( y_shift /= 0 ) THEN
	239
	240	IF ( bc_lr /= 'cyclic' .OR. bc_ns /= 'cyclic' ) THEN
	241	message_string = 'y_shift /= 0 is only allowed for cyclic ' // &
	242	'boundary conditions in both directions '
	243	CALL message( 'check_parameters', 'PA0467', 1, 2, 0, 6, 0 )
	244	ENDIF
	245	IF ( TRIM( psolver ) /= 'multigrid' .AND. &
	246	TRIM( psolver ) /= 'multigrid_noopt') &
	247	THEN
	248	message_string = 'y_shift /= 0 requires a multigrid pressure solver '
	249	CALL message( 'check_parameters', 'PA0468', 1, 2, 0, 6, 0 )
	250	ENDIF
	251
	252	CALL MPI_CART_COORDS( comm2d, pright, ndim, rcoord, ierr )
	253	CALL MPI_CART_COORDS( comm2d, pleft, ndim, lcoord, ierr )
	254
	255	!
	256	!-- If the x(y)-coordinate of the right (left) neighbor is smaller (greater)
	257	!-- than that of the calling process, then the calling process is located on
	258	!-- the right (left) boundary of the processor grid. In that case,
	259	!-- the y-coordinate of that neighbor is increased (decreased) by y_shift.
	260	!-- The rank of the process with that coordinate is then inquired and the
	261	!-- neighbor rank for MPI_SENDRECV, pright (pleft) is set to it.
	262	!-- In this way, the calling process receives a new right (left) neighbor
	263	!-- for all future MPI_SENDRECV calls. That neighbor has a y-coordinate
	264	!-- of y+(-)y_shift, where y is the original right (left) neighbor's
	265	!-- y-coordinate. The modulo-operation ensures that if the neighbor's
	266	!-- y-coordinate exceeds the grid-boundary, it will be relocated to
	267	!-- the opposite part of the grid cyclicly.
	268	IF ( rcoord(1) < pcoord(1) ) THEN
	269	rcoord(2) = MODULO( rcoord(2) + y_shift, pdims(2) )
	270	CALL MPI_CART_RANK( comm2d, rcoord, pright, ierr )
	271	ENDIF
	272
	273	IF ( lcoord(1) > pcoord(1) ) THEN
	274	lcoord(2) = MODULO( lcoord(2) - y_shift, pdims(2) )
	275	CALL MPI_CART_RANK( comm2d, lcoord, pleft, ierr )
	276	ENDIF
	277	ENDIF
	278	!
[2365]	279	!-- Vertical nesting: store four lists that identify partner ranks to exchange
	280	!-- data
	281	IF ( vnested ) CALL vnest_init_pegrid_rank
	282
	283	!
[1]	284	!-- Determine sub-topologies for transpositions
	285	!-- Transposition from z to x:
	286	remain_dims(1) = .TRUE.
	287	remain_dims(2) = .FALSE.
	288	CALL MPI_CART_SUB( comm2d, remain_dims, comm1dx, ierr )
	289	CALL MPI_COMM_RANK( comm1dx, myidx, ierr )
	290	!
	291	!-- Transposition from x to y
	292	remain_dims(1) = .FALSE.
	293	remain_dims(2) = .TRUE.
	294	CALL MPI_CART_SUB( comm2d, remain_dims, comm1dy, ierr )
	295	CALL MPI_COMM_RANK( comm1dy, myidy, ierr )
	296
	297
	298	!
[1003]	299	!-- Calculate array bounds along x-direction for every PE.
[3045]	300	ALLOCATE( nxlf(0:pdims(1)-1), nxrf(0:pdims(1)-1), nynf(0:pdims(2)-1), &
[1003]	301	nysf(0:pdims(2)-1) )
[1]	302
[1003]	303	IF ( MOD( nx+1 , pdims(1) ) /= 0 ) THEN
[3045]	304	WRITE( message_string, * ) 'x-direction: gridpoint number (',nx+1,') ', &
[4264]	305	'is not an& integral multiple of the number',&
	306	' of processors (',pdims(1),')'
[254]	307	CALL message( 'init_pegrid', 'PA0225', 1, 2, 0, 6, 0 )
[1]	308	ELSE
[1003]	309	nnx = ( nx + 1 ) / pdims(1)
[4264]	310	ENDIF
[1]	311
	312	!
	313	!-- Left and right array bounds, number of gridpoints
	314	DO i = 0, pdims(1)-1
	315	nxlf(i) = i * nnx
	316	nxrf(i) = ( i + 1 ) * nnx - 1
	317	ENDDO
	318
	319	!
	320	!-- Calculate array bounds in y-direction for every PE.
[1003]	321	IF ( MOD( ny+1 , pdims(2) ) /= 0 ) THEN
[274]	322	WRITE( message_string, * ) 'y-direction: gridpoint number (',ny+1,') ', &
[4264]	323	'is not an& integral multiple of the number',&
	324	' of processors (',pdims(2),')'
[254]	325	CALL message( 'init_pegrid', 'PA0227', 1, 2, 0, 6, 0 )
[1]	326	ELSE
[1003]	327	nny = ( ny + 1 ) / pdims(2)
[4264]	328	ENDIF
[1]	329
	330	!
	331	!-- South and north array bounds
	332	DO j = 0, pdims(2)-1
	333	nysf(j) = j * nny
	334	nynf(j) = ( j + 1 ) * nny - 1
	335	ENDDO
	336
	337	!
	338	!-- Local array bounds of the respective PEs
[1003]	339	nxl = nxlf(pcoord(1))
	340	nxr = nxrf(pcoord(1))
	341	nys = nysf(pcoord(2))
	342	nyn = nynf(pcoord(2))
	343	nzb = 0
	344	nzt = nz
	345	nnz = nz
[1]	346
	347	!
[707]	348	!-- Set switches to define if the PE is situated at the border of the virtual
	349	!-- processor grid
	350	IF ( nxl == 0 ) left_border_pe = .TRUE.
	351	IF ( nxr == nx ) right_border_pe = .TRUE.
	352	IF ( nys == 0 ) south_border_pe = .TRUE.
	353	IF ( nyn == ny ) north_border_pe = .TRUE.
	354
	355	!
[1]	356	!-- Calculate array bounds and gridpoint numbers for the transposed arrays
	357	!-- (needed in the pressure solver)
	358	!-- For the transposed arrays, cyclic boundaries as well as top and bottom
	359	!-- boundaries are omitted, because they are obstructive to the transposition
	360
	361	!
	362	!-- 1. transposition z --> x
[1001]	363	!-- This transposition is not neccessary in case of a 1d-decomposition along x
[2938]	364	IF ( psolver == 'poisfft' .OR. calculate_spectra ) THEN
[1304]	365
[1922]	366	IF ( pdims(2) /= 1 ) THEN
	367	IF ( MOD( nz , pdims(1) ) /= 0 ) THEN
[4264]	368	WRITE( message_string, * ) 'transposition z --> x:& ', &
	369	'nz=',nz,' is not an integral multiple ',&
	370	'of pdims(1)=',pdims(1)
[1922]	371	CALL message( 'init_pegrid', 'PA0230', 1, 2, 0, 6, 0 )
	372	ENDIF
[1]	373	ENDIF
[1922]	374
	375	nys_x = nys
	376	nyn_x = nyn
	377	nny_x = nny
	378	nnz_x = nz / pdims(1)
	379	nzb_x = 1 + myidx * nnz_x
	380	nzt_x = ( myidx + 1 ) * nnz_x
	381	sendrecvcount_zx = nnx * nny * nnz_x
	382
[1]	383	ENDIF
	384
[1922]	385
[4264]	386	IF ( psolver == 'poisfft' ) THEN
[1]	387	!
[1922]	388	!-- 2. transposition x --> y
	389	IF ( MOD( nx+1 , pdims(2) ) /= 0 ) THEN
[4264]	390	WRITE( message_string, * ) 'transposition x --> y:& ', &
	391	'nx+1=',nx+1,' is not an integral ', &
	392	'multiple of pdims(2)=',pdims(2)
[1922]	393	CALL message( 'init_pegrid', 'PA0231', 1, 2, 0, 6, 0 )
	394	ENDIF
[1]	395
[1922]	396	nnz_y = nnz_x
	397	nzb_y = nzb_x
	398	nzt_y = nzt_x
	399	nnx_y = (nx+1) / pdims(2)
	400	nxl_y = myidy * nnx_y
	401	nxr_y = ( myidy + 1 ) * nnx_y - 1
	402	sendrecvcount_xy = nnx_y * nny_x * nnz_y
[1]	403	!
[4264]	404	!-- 3. transposition y --> z
[1922]	405	!-- (ELSE: x --> y in case of 1D-decomposition along x)
	406	nxl_z = nxl_y
	407	nxr_z = nxr_y
	408	nny_z = (ny+1) / pdims(1)
	409	nys_z = myidx * nny_z
	410	nyn_z = ( myidx + 1 ) * nny_z - 1
	411	sendrecvcount_yz = nnx_y * nny_z * nnz_y
[1304]	412
[1922]	413	IF ( pdims(2) /= 1 ) THEN
[1]	414	!
[1922]	415	!-- y --> z
	416	!-- This transposition is not neccessary in case of a 1d-decomposition
	417	!-- along x, except that the uptream-spline method is switched on
	418	IF ( MOD( ny+1 , pdims(1) ) /= 0 ) THEN
[4264]	419	WRITE( message_string, * ) 'transposition y --> z:& ', &
	420	'ny+1=',ny+1,' is not an integral ', &
	421	'multiple of pdims(1)=',pdims(1)
[1922]	422	CALL message( 'init_pegrid', 'PA0232', 1, 2, 0, 6, 0 )
	423	ENDIF
[1]	424
[1922]	425	ELSE
[1]	426	!
[1922]	427	!-- x --> y
	428	!-- This condition must be fulfilled for a 1D-decomposition along x
	429	IF ( MOD( ny+1 , pdims(1) ) /= 0 ) THEN
[4264]	430	WRITE( message_string, * ) 'transposition x --> y:& ', &
	431	'ny+1=',ny+1,' is not an integral ', &
	432	'multiple of pdims(1)=',pdims(1)
[1922]	433	CALL message( 'init_pegrid', 'PA0233', 1, 2, 0, 6, 0 )
	434	ENDIF
	435
[1]	436	ENDIF
	437
	438	ENDIF
	439
	440	!
	441	!-- Indices for direct transpositions z --> y (used for calculating spectra)
[1922]	442	IF ( calculate_spectra ) THEN
[1003]	443	IF ( MOD( nz, pdims(2) ) /= 0 ) THEN
[1922]	444	WRITE( message_string, * ) 'direct transposition z --> y (needed ', &
[4264]	445	'for spectra):& nz=',nz,' is not an ', &
	446	'integral multiple of pdims(2)=',pdims(2)
[254]	447	CALL message( 'init_pegrid', 'PA0234', 1, 2, 0, 6, 0 )
[1]	448	ELSE
[1003]	449	nxl_yd = nxl
	450	nxr_yd = nxr
	451	nzb_yd = 1 + myidy * ( nz / pdims(2) )
	452	nzt_yd = ( myidy + 1 ) * ( nz / pdims(2) )
	453	sendrecvcount_zyd = nnx * nny * ( nz / pdims(2) )
[1]	454	ENDIF
	455	ENDIF
	456
[1922]	457	IF ( psolver == 'poisfft' .OR. calculate_spectra ) THEN
[1]	458	!
[1922]	459	!-- Indices for direct transpositions y --> x
	460	!-- (they are only possible in case of a 1d-decomposition along x)
	461	IF ( pdims(2) == 1 ) THEN
	462	nny_x = nny / pdims(1)
	463	nys_x = myid * nny_x
	464	nyn_x = ( myid + 1 ) * nny_x - 1
	465	nzb_x = 1
	466	nzt_x = nz
	467	sendrecvcount_xy = nnx * nny_x * nz
	468	ENDIF
	469
[1]	470	ENDIF
	471
[1922]	472	IF ( psolver == 'poisfft' ) THEN
[1]	473	!
[1922]	474	!-- Indices for direct transpositions x --> y
	475	!-- (they are only possible in case of a 1d-decomposition along y)
	476	IF ( pdims(1) == 1 ) THEN
	477	nnx_y = nnx / pdims(2)
	478	nxl_y = myid * nnx_y
	479	nxr_y = ( myid + 1 ) * nnx_y - 1
	480	nzb_y = 1
	481	nzt_y = nz
	482	sendrecvcount_xy = nnx_y * nny * nz
	483	ENDIF
	484
[1]	485	ENDIF
	486
	487	!
	488	!-- Arrays for storing the array bounds are needed any more
	489	DEALLOCATE( nxlf , nxrf , nynf , nysf )
	490
[807]	491
[145]	492	!
	493	!-- Collect index bounds from other PEs (to be written to restart file later)
	494	ALLOCATE( hor_index_bounds(4,0:numprocs-1) )
	495
	496	IF ( myid == 0 ) THEN
	497
	498	hor_index_bounds(1,0) = nxl
	499	hor_index_bounds(2,0) = nxr
	500	hor_index_bounds(3,0) = nys
	501	hor_index_bounds(4,0) = nyn
	502
	503	!
	504	!-- Receive data from all other PEs
	505	DO i = 1, numprocs-1
	506	CALL MPI_RECV( ibuf, 4, MPI_INTEGER, i, MPI_ANY_TAG, comm2d, status, &
	507	ierr )
	508	hor_index_bounds(:,i) = ibuf(1:4)
	509	ENDDO
	510
	511	ELSE
	512	!
	513	!-- Send index bounds to PE0
	514	ibuf(1) = nxl
	515	ibuf(2) = nxr
	516	ibuf(3) = nys
	517	ibuf(4) = nyn
	518	CALL MPI_SEND( ibuf, 4, MPI_INTEGER, 0, myid, comm2d, ierr )
	519
	520	ENDIF
	521
[807]	522
[1]	523	#if defined( __print )
	524	!
	525	!-- Control output
	526	IF ( myid == 0 ) THEN
	527	PRINT, ' processor topology *'
	528	PRINT*, ' '
	529	PRINT*, 'myid pcoord left right south north idx idy nxl: nxr',&
	530	&' nys: nyn'
	531	PRINT*, '------------------------------------------------------------',&
	532	&'-----------'
	533	WRITE (*,1000) 0, pcoord(1), pcoord(2), pleft, pright, psouth, pnorth, &
	534	myidx, myidy, nxl, nxr, nys, nyn
	535	1000 FORMAT (I4,2X,'(',I3,',',I3,')',3X,I4,2X,I4,3X,I4,2X,I4,2X,I3,1X,I3, &
	536	2(2X,I4,':',I4))
	537
	538	!
[108]	539	!-- Receive data from the other PEs
[1]	540	DO i = 1,numprocs-1
	541	CALL MPI_RECV( ibuf, 12, MPI_INTEGER, i, MPI_ANY_TAG, comm2d, status, &
	542	ierr )
	543	WRITE (*,1000) i, ( ibuf(j) , j = 1,12 )
	544	ENDDO
	545	ELSE
	546
	547	!
	548	!-- Send data to PE0
	549	ibuf(1) = pcoord(1); ibuf(2) = pcoord(2); ibuf(3) = pleft
	550	ibuf(4) = pright; ibuf(5) = psouth; ibuf(6) = pnorth; ibuf(7) = myidx
	551	ibuf(8) = myidy; ibuf(9) = nxl; ibuf(10) = nxr; ibuf(11) = nys
	552	ibuf(12) = nyn
	553	CALL MPI_SEND( ibuf, 12, MPI_INTEGER, 0, myid, comm2d, ierr )
	554	ENDIF
	555	#endif
	556
[667]	557	!
[709]	558	!-- Determine the number of ghost point layers
[3347]	559	IF ( scalar_advec == 'ws-scheme' .OR. &
[3999]	560	momentum_advec == 'ws-scheme' .OR. nested_run ) THEN
[667]	561	nbgp = 3
	562	ELSE
	563	nbgp = 1
[709]	564	ENDIF
[667]	565
[102]	566	!
[4241]	567	!-- Check that the number of computational grid points is not smaller than the number of
	568	!-- ghost points.
	569	IF ( nnx < nbgp ) THEN
	570	WRITE( message_string, * ) 'number of subdomain grid points along x (', nnx, ') is smaller',&
	571	'than the number of ghost points (', nbgp, ')'
	572	CALL message( 'init_pegrid', 'PA0682', 1, 2, 0, 6, 0 )
	573	ENDIF
	574	IF ( nny < nbgp ) THEN
	575	WRITE( message_string, * ) 'number of subdomain grid points along y (', nny, ') is smaller',&
	576	'than the number of ghost points (', nbgp, ')'
	577	CALL message( 'init_pegrid', 'PA0683', 1, 2, 0, 6, 0 )
	578	ENDIF
	579
	580	!
[709]	581	!-- Create a new MPI derived datatype for the exchange of surface (xy) data,
	582	!-- which is needed for coupled atmosphere-ocean runs.
	583	!-- First, calculate number of grid points of an xy-plane.
[667]	584	ngp_xy = ( nxr - nxl + 1 + 2 * nbgp ) * ( nyn - nys + 1 + 2 * nbgp )
[102]	585	CALL MPI_TYPE_VECTOR( ngp_xy, 1, nzt-nzb+2, MPI_REAL, type_xy, ierr )
	586	CALL MPI_TYPE_COMMIT( type_xy, ierr )
[667]	587
[2365]	588	IF ( TRIM( coupling_mode ) /= 'uncoupled' .AND. .NOT. vnested ) THEN
[667]	589
	590	!
	591	!-- Pass the number of grid points of the atmosphere model to
	592	!-- the ocean model and vice versa
	593	IF ( coupling_mode == 'atmosphere_to_ocean' ) THEN
	594
	595	nx_a = nx
	596	ny_a = ny
	597
[709]	598	IF ( myid == 0 ) THEN
	599
	600	CALL MPI_SEND( nx_a, 1, MPI_INTEGER, numprocs, 1, comm_inter, &
	601	ierr )
	602	CALL MPI_SEND( ny_a, 1, MPI_INTEGER, numprocs, 2, comm_inter, &
	603	ierr )
	604	CALL MPI_SEND( pdims, 2, MPI_INTEGER, numprocs, 3, comm_inter, &
	605	ierr )
	606	CALL MPI_RECV( nx_o, 1, MPI_INTEGER, numprocs, 4, comm_inter, &
	607	status, ierr )
	608	CALL MPI_RECV( ny_o, 1, MPI_INTEGER, numprocs, 5, comm_inter, &
	609	status, ierr )
	610	CALL MPI_RECV( pdims_remote, 2, MPI_INTEGER, numprocs, 6, &
[667]	611	comm_inter, status, ierr )
	612	ENDIF
	613
[709]	614	CALL MPI_BCAST( nx_o, 1, MPI_INTEGER, 0, comm2d, ierr )
	615	CALL MPI_BCAST( ny_o, 1, MPI_INTEGER, 0, comm2d, ierr )
	616	CALL MPI_BCAST( pdims_remote, 2, MPI_INTEGER, 0, comm2d, ierr )
[667]	617
	618	ELSEIF ( coupling_mode == 'ocean_to_atmosphere' ) THEN
	619
	620	nx_o = nx
	621	ny_o = ny
	622
	623	IF ( myid == 0 ) THEN
[709]	624
	625	CALL MPI_RECV( nx_a, 1, MPI_INTEGER, 0, 1, comm_inter, status, &
	626	ierr )
	627	CALL MPI_RECV( ny_a, 1, MPI_INTEGER, 0, 2, comm_inter, status, &
	628	ierr )
	629	CALL MPI_RECV( pdims_remote, 2, MPI_INTEGER, 0, 3, comm_inter, &
	630	status, ierr )
	631	CALL MPI_SEND( nx_o, 1, MPI_INTEGER, 0, 4, comm_inter, ierr )
	632	CALL MPI_SEND( ny_o, 1, MPI_INTEGER, 0, 5, comm_inter, ierr )
	633	CALL MPI_SEND( pdims, 2, MPI_INTEGER, 0, 6, comm_inter, ierr )
[667]	634	ENDIF
	635
	636	CALL MPI_BCAST( nx_a, 1, MPI_INTEGER, 0, comm2d, ierr)
	637	CALL MPI_BCAST( ny_a, 1, MPI_INTEGER, 0, comm2d, ierr)
	638	CALL MPI_BCAST( pdims_remote, 2, MPI_INTEGER, 0, comm2d, ierr)
	639
	640	ENDIF
	641
[709]	642	ngp_a = ( nx_a+1 + 2 * nbgp ) * ( ny_a+1 + 2 * nbgp )
	643	ngp_o = ( nx_o+1 + 2 * nbgp ) * ( ny_o+1 + 2 * nbgp )
[667]	644
	645	!
[709]	646	!-- Determine if the horizontal grid and the number of PEs in ocean and
	647	!-- atmosphere is same or not
	648	IF ( nx_o == nx_a .AND. ny_o == ny_a .AND. &
[667]	649	pdims(1) == pdims_remote(1) .AND. pdims(2) == pdims_remote(2) ) &
	650	THEN
	651	coupling_topology = 0
	652	ELSE
	653	coupling_topology = 1
	654	ENDIF
	655
	656	!
	657	!-- Determine the target PEs for the exchange between ocean and
	658	!-- atmosphere (comm2d)
[709]	659	IF ( coupling_topology == 0 ) THEN
	660	!
	661	!-- In case of identical topologies, every atmosphere PE has exactly one
	662	!-- ocean PE counterpart and vice versa
	663	IF ( TRIM( coupling_mode ) == 'atmosphere_to_ocean' ) THEN
[667]	664	target_id = myid + numprocs
	665	ELSE
	666	target_id = myid
	667	ENDIF
	668
	669	ELSE
	670	!
	671	!-- In case of nonequivalent topology in ocean and atmosphere only for
	672	!-- PE0 in ocean and PE0 in atmosphere a target_id is needed, since
[709]	673	!-- data echxchange between ocean and atmosphere will be done only
	674	!-- between these PEs.
	675	IF ( myid == 0 ) THEN
	676
	677	IF ( TRIM( coupling_mode ) == 'atmosphere_to_ocean' ) THEN
[667]	678	target_id = numprocs
	679	ELSE
	680	target_id = 0
	681	ENDIF
[709]	682
[667]	683	ENDIF
[709]	684
[667]	685	ENDIF
	686
	687	ENDIF
	688
[2365]	689	!
	690	!-- Store partner grid point co-ordinates as lists.
	691	!-- Create custom MPI vector datatypes for contiguous data transfer
	692	IF ( vnested ) CALL vnest_init_pegrid_domain
[667]	693
[1]	694	#else
	695
	696	!
	697	!-- Array bounds when running on a single PE (respectively a non-parallel
	698	!-- machine)
[1003]	699	nxl = 0
	700	nxr = nx
	701	nnx = nxr - nxl + 1
	702	nys = 0
	703	nyn = ny
	704	nny = nyn - nys + 1
	705	nzb = 0
	706	nzt = nz
	707	nnz = nz
[1]	708
[145]	709	ALLOCATE( hor_index_bounds(4,0:0) )
	710	hor_index_bounds(1,0) = nxl
	711	hor_index_bounds(2,0) = nxr
	712	hor_index_bounds(3,0) = nys
	713	hor_index_bounds(4,0) = nyn
	714
[1]	715	!
	716	!-- Array bounds for the pressure solver (in the parallel code, these bounds
	717	!-- are the ones for the transposed arrays)
[1003]	718	nys_x = nys
	719	nyn_x = nyn
	720	nzb_x = nzb + 1
	721	nzt_x = nzt
[1]	722
[1003]	723	nxl_y = nxl
	724	nxr_y = nxr
	725	nzb_y = nzb + 1
	726	nzt_y = nzt
[1]	727
[1003]	728	nxl_z = nxl
	729	nxr_z = nxr
	730	nys_z = nys
	731	nyn_z = nyn
[1]	732
	733	#endif
	734
	735	!
	736	!-- Calculate number of grid levels necessary for the multigrid poisson solver
	737	!-- as well as the gridpoint indices on each level
[1575]	738	IF ( psolver(1:9) == 'multigrid' ) THEN
[1]	739
	740	!
	741	!-- First calculate number of possible grid levels for the subdomains
	742	mg_levels_x = 1
	743	mg_levels_y = 1
	744	mg_levels_z = 1
	745
	746	i = nnx
	747	DO WHILE ( MOD( i, 2 ) == 0 .AND. i /= 2 )
	748	i = i / 2
	749	mg_levels_x = mg_levels_x + 1
	750	ENDDO
	751
	752	j = nny
	753	DO WHILE ( MOD( j, 2 ) == 0 .AND. j /= 2 )
	754	j = j / 2
	755	mg_levels_y = mg_levels_y + 1
	756	ENDDO
	757
[181]	758	k = nz ! do not use nnz because it might be > nz due to transposition
	759	! requirements
[1]	760	DO WHILE ( MOD( k, 2 ) == 0 .AND. k /= 2 )
	761	k = k / 2
	762	mg_levels_z = mg_levels_z + 1
	763	ENDDO
[2197]	764	!
	765	!-- The optimized MG-solver does not allow odd values for nz at the coarsest
	766	!-- grid level
	767	IF ( TRIM( psolver ) /= 'multigrid_noopt' ) THEN
	768	IF ( MOD( k, 2 ) /= 0 ) mg_levels_z = mg_levels_z - 1
[3057]	769	!
	770	!-- An odd value of nz does not work. The finest level must have an even
	771	!-- value.
	772	IF ( mg_levels_z == 0 ) THEN
	773	message_string = 'optimized multigrid method requires nz to be even'
[3058]	774	CALL message( 'init_pegrid', 'PA0495', 1, 2, 0, 6, 0 )
[3057]	775	ENDIF
[2197]	776	ENDIF
[1]	777
	778	maximum_grid_level = MIN( mg_levels_x, mg_levels_y, mg_levels_z )
[3890]	779	!
	780	!-- Check if subdomain sizes prevents any coarsening.
	781	!-- This case, the maximum number of grid levels is 1, i.e. effectively
	782	!-- a Gauss-Seidel scheme is applied rather than a multigrid approach.
[3897]	783	!-- Give a warning in this case.
	784	IF ( maximum_grid_level == 1 .AND. mg_switch_to_pe0_level == -1 ) THEN
	785	message_string = 'No grid coarsening possible, multigrid ' // &
	786	'approach effectively reduces to a Gauss-Seidel ' //&
	787	'scheme.'
	788
	789	CALL message( 'poismg', 'PA0648', 0, 1, 0, 6, 0 )
[3890]	790	ENDIF
[1]	791
	792	!
	793	!-- Find out, if the total domain allows more levels. These additional
[709]	794	!-- levels are identically processed on all PEs.
[197]	795	IF ( numprocs > 1 .AND. mg_switch_to_pe0_level /= -1 ) THEN
[709]	796
[1]	797	IF ( mg_levels_z > MIN( mg_levels_x, mg_levels_y ) ) THEN
[709]	798
[1]	799	mg_switch_to_pe0_level_l = maximum_grid_level
	800
	801	mg_levels_x = 1
	802	mg_levels_y = 1
	803
	804	i = nx+1
	805	DO WHILE ( MOD( i, 2 ) == 0 .AND. i /= 2 )
	806	i = i / 2
	807	mg_levels_x = mg_levels_x + 1
	808	ENDDO
	809
	810	j = ny+1
	811	DO WHILE ( MOD( j, 2 ) == 0 .AND. j /= 2 )
	812	j = j / 2
	813	mg_levels_y = mg_levels_y + 1
	814	ENDDO
	815
	816	maximum_grid_level_l = MIN( mg_levels_x, mg_levels_y, mg_levels_z )
	817
	818	IF ( maximum_grid_level_l > mg_switch_to_pe0_level_l ) THEN
	819	mg_switch_to_pe0_level_l = maximum_grid_level_l - &
	820	mg_switch_to_pe0_level_l + 1
	821	ELSE
	822	mg_switch_to_pe0_level_l = 0
	823	ENDIF
[709]	824
[1]	825	ELSE
[3057]	826
[1]	827	mg_switch_to_pe0_level_l = 0
	828	maximum_grid_level_l = maximum_grid_level
[709]	829
[1]	830	ENDIF
	831
	832	!
	833	!-- Use switch level calculated above only if it is not pre-defined
	834	!-- by user
	835	IF ( mg_switch_to_pe0_level == 0 ) THEN
	836	IF ( mg_switch_to_pe0_level_l /= 0 ) THEN
	837	mg_switch_to_pe0_level = mg_switch_to_pe0_level_l
	838	maximum_grid_level = maximum_grid_level_l
	839	ENDIF
	840
	841	ELSE
	842	!
	843	!-- Check pre-defined value and reset to default, if neccessary
[3045]	844	IF ( mg_switch_to_pe0_level < mg_switch_to_pe0_level_l .OR. &
[1]	845	mg_switch_to_pe0_level >= maximum_grid_level_l ) THEN
[3045]	846	message_string = 'mg_switch_to_pe0_level ' // &
[2271]	847	'out of range and reset to 0'
[254]	848	CALL message( 'init_pegrid', 'PA0235', 0, 1, 0, 6, 0 )
[1]	849	mg_switch_to_pe0_level = 0
	850	ELSE
	851	!
	852	!-- Use the largest number of possible levels anyway and recalculate
	853	!-- the switch level to this largest number of possible values
	854	maximum_grid_level = maximum_grid_level_l
	855
	856	ENDIF
[709]	857
[1]	858	ENDIF
	859
	860	ENDIF
	861
[1056]	862	ALLOCATE( grid_level_count(maximum_grid_level), &
	863	nxl_mg(0:maximum_grid_level), nxr_mg(0:maximum_grid_level), &
	864	nyn_mg(0:maximum_grid_level), nys_mg(0:maximum_grid_level), &
	865	nzt_mg(0:maximum_grid_level) )
[1]	866
	867	grid_level_count = 0
[1056]	868	!
	869	!-- Index zero required as dummy due to definition of arrays f2 and p2 in
	870	!-- recursive subroutine next_mg_level
	871	nxl_mg(0) = 0; nxr_mg(0) = 0; nyn_mg(0) = 0; nys_mg(0) = 0; nzt_mg(0) = 0
[778]	872
[1]	873	nxl_l = nxl; nxr_l = nxr; nys_l = nys; nyn_l = nyn; nzt_l = nzt
	874
	875	DO i = maximum_grid_level, 1 , -1
	876
	877	IF ( i == mg_switch_to_pe0_level ) THEN
[1804]	878	#if defined( __parallel )
[1]	879	!
	880	!-- Save the grid size of the subdomain at the switch level, because
	881	!-- it is needed in poismg.
	882	ind(1) = nxl_l; ind(2) = nxr_l
	883	ind(3) = nys_l; ind(4) = nyn_l
	884	ind(5) = nzt_l
	885	ALLOCATE( ind_all(5*numprocs), mg_loc_ind(5,0:numprocs-1) )
	886	CALL MPI_ALLGATHER( ind, 5, MPI_INTEGER, ind_all, 5, &
	887	MPI_INTEGER, comm2d, ierr )
	888	DO j = 0, numprocs-1
	889	DO k = 1, 5
	890	mg_loc_ind(k,j) = ind_all(k+j*5)
	891	ENDDO
	892	ENDDO
	893	DEALLOCATE( ind_all )
	894	!
[709]	895	!-- Calculate the grid size of the total domain
[1]	896	nxr_l = ( nxr_l-nxl_l+1 ) * pdims(1) - 1
	897	nxl_l = 0
	898	nyn_l = ( nyn_l-nys_l+1 ) * pdims(2) - 1
	899	nys_l = 0
	900	!
	901	!-- The size of this gathered array must not be larger than the
	902	!-- array tend, which is used in the multigrid scheme as a temporary
[778]	903	!-- array. Therefore the subdomain size of an PE is calculated and
	904	!-- the size of the gathered grid. These values are used in
	905	!-- routines pres and poismg
	906	subdomain_size = ( nxr - nxl + 2 * nbgp + 1 ) * &
	907	( nyn - nys + 2 * nbgp + 1 ) * ( nzt - nzb + 2 )
[3045]	908	gathered_size = ( nxr_l - nxl_l + 3 ) * ( nyn_l - nys_l + 3 ) * &
[1]	909	( nzt_l - nzb + 2 )
	910
[1804]	911	#else
[3045]	912	message_string = 'multigrid gather/scatter impossible ' // &
[1]	913	'in non parallel mode'
[254]	914	CALL message( 'init_pegrid', 'PA0237', 1, 2, 0, 6, 0 )
[1]	915	#endif
	916	ENDIF
	917
	918	nxl_mg(i) = nxl_l
	919	nxr_mg(i) = nxr_l
	920	nys_mg(i) = nys_l
	921	nyn_mg(i) = nyn_l
	922	nzt_mg(i) = nzt_l
	923
	924	nxl_l = nxl_l / 2
	925	nxr_l = nxr_l / 2
	926	nys_l = nys_l / 2
	927	nyn_l = nyn_l / 2
	928	nzt_l = nzt_l / 2
[3890]	929
[1]	930	ENDDO
	931
[780]	932	!
[3045]	933	!-- Temporary problem: Currently calculation of maxerror in routine poismg crashes
[780]	934	!-- if grid data are collected on PE0 already on the finest grid level.
	935	!-- To be solved later.
	936	IF ( maximum_grid_level == mg_switch_to_pe0_level ) THEN
	937	message_string = 'grid coarsening on subdomain level cannot be performed'
	938	CALL message( 'poismg', 'PA0236', 1, 2, 0, 6, 0 )
	939	ENDIF
	940
[1]	941	ELSE
	942
[667]	943	maximum_grid_level = 0
[1]	944
	945	ENDIF
	946
[722]	947	!
	948	!-- Default level 0 tells exchange_horiz that all ghost planes have to be
	949	!-- exchanged. grid_level is adjusted in poismg, where only one ghost plane
	950	!-- is required.
	951	grid_level = 0
[1]	952
[1804]	953	#if defined( __parallel )
[1]	954	!
	955	!-- Gridpoint number for the exchange of ghost points (y-line for 2D-arrays)
[667]	956	ngp_y = nyn - nys + 1 + 2 * nbgp
[1]	957
	958	!
[709]	959	!-- Define new MPI derived datatypes for the exchange of ghost points in
	960	!-- x- and y-direction for 2D-arrays (line)
[1968]	961	CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp, ngp_y, MPI_REAL, type_x, &
[709]	962	ierr )
[1]	963	CALL MPI_TYPE_COMMIT( type_x, ierr )
	964
[667]	965	CALL MPI_TYPE_VECTOR( nbgp, ngp_y, ngp_y, MPI_REAL, type_y, ierr )
	966	CALL MPI_TYPE_COMMIT( type_y, ierr )
[1968]	967	!
	968	!-- Define new MPI derived datatypes for the exchange of ghost points in
[3542]	969	!-- x- and y-direction for 2D-INTEGER arrays (line) - on normal grid.
	970	!-- Define types for 32-bit and 8-bit Integer. The 8-bit Integer are only
	971	!-- required on normal grid, while 32-bit Integer may be also required on
	972	!-- coarser grid level in case of multigrid solver.
	973	!
	974	!-- 8-bit Integer
	975	CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp, ngp_y, MPI_BYTE, &
	976	type_x_byte, ierr )
	977	CALL MPI_TYPE_COMMIT( type_x_byte, ierr )
	978
	979	CALL MPI_TYPE_VECTOR( nbgp, ngp_y, ngp_y, MPI_BYTE, &
	980	type_y_byte, ierr )
	981	CALL MPI_TYPE_COMMIT( type_y_byte, ierr )
	982	!
	983	!-- 32-bit Integer
[1968]	984	ALLOCATE( type_x_int(0:maximum_grid_level), &
	985	type_y_int(0:maximum_grid_level) )
[3542]	986
[1968]	987	CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp, ngp_y, MPI_INTEGER, &
	988	type_x_int(0), ierr )
	989	CALL MPI_TYPE_COMMIT( type_x_int(0), ierr )
[667]	990
[1968]	991	CALL MPI_TYPE_VECTOR( nbgp, ngp_y, ngp_y, MPI_INTEGER, type_y_int(0), ierr )
	992	CALL MPI_TYPE_COMMIT( type_y_int(0), ierr )
[1]	993	!
	994	!-- Calculate gridpoint numbers for the exchange of ghost points along x
	995	!-- (yz-plane for 3D-arrays) and define MPI derived data type(s) for the
	996	!-- exchange of ghost points in y-direction (xz-plane).
	997	!-- Do these calculations for the model grid and (if necessary) also
	998	!-- for the coarser grid levels used in the multigrid method
[2696]	999	ALLOCATE ( ngp_xz(0:maximum_grid_level), &
	1000	ngp_xz_int(0:maximum_grid_level), &
	1001	ngp_yz(0:maximum_grid_level), &
	1002	ngp_yz_int(0:maximum_grid_level), &
	1003	type_xz(0:maximum_grid_level), &
	1004	type_xz_int(0:maximum_grid_level), &
	1005	type_yz(0:maximum_grid_level), &
	1006	type_yz_int(0:maximum_grid_level) )
[1]	1007
	1008	nxl_l = nxl; nxr_l = nxr; nys_l = nys; nyn_l = nyn; nzb_l = nzb; nzt_l = nzt
[709]	1009
[667]	1010	!
	1011	!-- Discern between the model grid, which needs nbgp ghost points and
	1012	!-- grid levels for the multigrid scheme. In the latter case only one
	1013	!-- ghost point is necessary.
[709]	1014	!-- First definition of MPI-datatypes for exchange of ghost layers on normal
[667]	1015	!-- grid. The following loop is needed for data exchange in poismg.f90.
	1016	!
	1017	!-- Determine number of grid points of yz-layer for exchange
	1018	ngp_yz(0) = (nzt - nzb + 2) * (nyn - nys + 1 + 2 * nbgp)
[709]	1019
[667]	1020	!
[709]	1021	!-- Define an MPI-datatype for the exchange of left/right boundaries.
	1022	!-- Although data are contiguous in physical memory (which does not
	1023	!-- necessarily require an MPI-derived datatype), the data exchange between
	1024	!-- left and right PE's using the MPI-derived type is 10% faster than without.
[667]	1025	CALL MPI_TYPE_VECTOR( nxr-nxl+1+2nbgp, nbgp(nzt-nzb+2), ngp_yz(0), &
[709]	1026	MPI_REAL, type_xz(0), ierr )
[667]	1027	CALL MPI_TYPE_COMMIT( type_xz(0), ierr )
[1]	1028
[709]	1029	CALL MPI_TYPE_VECTOR( nbgp, ngp_yz(0), ngp_yz(0), MPI_REAL, type_yz(0), &
	1030	ierr )
[667]	1031	CALL MPI_TYPE_COMMIT( type_yz(0), ierr )
[709]	1032
[667]	1033	!
[2696]	1034	!-- Define data types for exchange of 3D Integer arrays.
	1035	ngp_yz_int(0) = (nzt - nzb + 2) * (nyn - nys + 1 + 2 * nbgp)
	1036
	1037	CALL MPI_TYPE_VECTOR( nxr-nxl+1+2nbgp, nbgp(nzt-nzb+2), ngp_yz_int(0), &
	1038	MPI_INTEGER, type_xz_int(0), ierr )
	1039	CALL MPI_TYPE_COMMIT( type_xz_int(0), ierr )
	1040
	1041	CALL MPI_TYPE_VECTOR( nbgp, ngp_yz_int(0), ngp_yz_int(0), MPI_INTEGER, &
	1042	type_yz_int(0), ierr )
	1043	CALL MPI_TYPE_COMMIT( type_yz_int(0), ierr )
	1044
	1045	!
[709]	1046	!-- Definition of MPI-datatypes for multigrid method (coarser level grids)
[1575]	1047	IF ( psolver(1:9) == 'multigrid' ) THEN
[667]	1048	!
[709]	1049	!-- Definition of MPI-datatyoe as above, but only 1 ghost level is used
	1050	DO i = maximum_grid_level, 1 , -1
[1968]	1051	!
[2696]	1052	!-- For 3D-exchange on different multigrid level, one ghost point for
	1053	!-- REAL arrays, two ghost points for INTEGER arrays
[1575]	1054	ngp_xz(i) = (nzt_l - nzb_l + 2) * (nxr_l - nxl_l + 3)
[667]	1055	ngp_yz(i) = (nzt_l - nzb_l + 2) * (nyn_l - nys_l + 3)
	1056
[2696]	1057	ngp_xz_int(i) = (nzt_l - nzb_l + 2) * (nxr_l - nxl_l + 3)
	1058	ngp_yz_int(i) = (nzt_l - nzb_l + 2) * (nyn_l - nys_l + 3)
	1059	!
	1060	!-- MPI data type for REAL arrays, for xz-layers
	1061	CALL MPI_TYPE_VECTOR( nxr_l-nxl_l+3, nzt_l-nzb_l+2, ngp_yz(i), &
[709]	1062	MPI_REAL, type_xz(i), ierr )
[667]	1063	CALL MPI_TYPE_COMMIT( type_xz(i), ierr )
[1]	1064
[2696]	1065	!
	1066	!-- MPI data type for INTEGER arrays, for xz-layers
	1067	CALL MPI_TYPE_VECTOR( nxr_l-nxl_l+3, nzt_l-nzb_l+2, ngp_yz_int(i), &
	1068	MPI_INTEGER, type_xz_int(i), ierr )
	1069	CALL MPI_TYPE_COMMIT( type_xz_int(i), ierr )
	1070
	1071	!
	1072	!-- MPI data type for REAL arrays, for yz-layers
[709]	1073	CALL MPI_TYPE_VECTOR( 1, ngp_yz(i), ngp_yz(i), MPI_REAL, type_yz(i), &
	1074	ierr )
[667]	1075	CALL MPI_TYPE_COMMIT( type_yz(i), ierr )
[2696]	1076	!
	1077	!-- MPI data type for INTEGER arrays, for yz-layers
	1078	CALL MPI_TYPE_VECTOR( 1, ngp_yz_int(i), ngp_yz_int(i), MPI_INTEGER, &
	1079	type_yz_int(i), ierr )
	1080	CALL MPI_TYPE_COMMIT( type_yz_int(i), ierr )
[667]	1081
[1968]	1082
	1083	!-- For 2D-exchange of INTEGER arrays on coarser grid level, where 2 ghost
[3542]	1084	!-- points need to be exchanged. Only required for 32-bit Integer arrays.
	1085	CALL MPI_TYPE_VECTOR( nxr_l-nxl_l+5, 2, nyn_l-nys_l+5, MPI_INTEGER, &
[1968]	1086	type_x_int(i), ierr )
	1087	CALL MPI_TYPE_COMMIT( type_x_int(i), ierr )
	1088
	1089
[3542]	1090	CALL MPI_TYPE_VECTOR( 2, nyn_l-nys_l+5, nyn_l-nys_l+5, MPI_INTEGER, &
[1968]	1091	type_y_int(i), ierr )
	1092	CALL MPI_TYPE_COMMIT( type_y_int(i), ierr )
	1093
[667]	1094	nxl_l = nxl_l / 2
	1095	nxr_l = nxr_l / 2
	1096	nys_l = nys_l / 2
	1097	nyn_l = nyn_l / 2
	1098	nzt_l = nzt_l / 2
[709]	1099
[667]	1100	ENDDO
[709]	1101
	1102	ENDIF
[1677]	1103
[1]	1104	#endif
	1105
[1804]	1106	#if defined( __parallel )
[1]	1107	!
[1933]	1108	!-- Setting of flags for inflow/outflow/nesting conditions.
[1964]	1109	IF ( pleft == MPI_PROC_NULL ) THEN
[3182]	1110	IF ( bc_lr == 'dirichlet/radiation' .OR. bc_lr == 'nested' .OR. &
	1111	bc_lr == 'nesting_offline' ) THEN
	1112	bc_dirichlet_l = .TRUE.
[1964]	1113	ELSEIF ( bc_lr == 'radiation/dirichlet' ) THEN
[3182]	1114	bc_radiation_l = .TRUE.
[1]	1115	ENDIF
	1116	ENDIF
[1933]	1117
[1964]	1118	IF ( pright == MPI_PROC_NULL ) THEN
	1119	IF ( bc_lr == 'dirichlet/radiation' ) THEN
[3182]	1120	bc_radiation_r = .TRUE.
	1121	ELSEIF ( bc_lr == 'radiation/dirichlet' .OR. bc_lr == 'nested' .OR. &
	1122	bc_lr == 'nesting_offline' ) THEN
	1123	bc_dirichlet_r = .TRUE.
[1]	1124	ENDIF
	1125	ENDIF
	1126
[1964]	1127	IF ( psouth == MPI_PROC_NULL ) THEN
	1128	IF ( bc_ns == 'dirichlet/radiation' ) THEN
[3182]	1129	bc_radiation_s = .TRUE.
	1130	ELSEIF ( bc_ns == 'radiation/dirichlet' .OR. bc_ns == 'nested' .OR. &
	1131	bc_ns == 'nesting_offline' ) THEN
	1132	bc_dirichlet_s = .TRUE.
[1]	1133	ENDIF
	1134	ENDIF
	1135
[1964]	1136	IF ( pnorth == MPI_PROC_NULL ) THEN
[3182]	1137	IF ( bc_ns == 'dirichlet/radiation' .OR. bc_ns == 'nested' .OR. &
	1138	bc_ns == 'nesting_offline' ) THEN
	1139	bc_dirichlet_n = .TRUE.
[1964]	1140	ELSEIF ( bc_ns == 'radiation/dirichlet' ) THEN
[3182]	1141	bc_radiation_n = .TRUE.
[1]	1142	ENDIF
	1143	ENDIF
[2938]	1144	!
	1145	!-- In case of synthetic turbulence geneartor determine ids.
	1146	!-- Please note, if no forcing or nesting is applied, the generator is applied
	1147	!-- only at the left lateral boundary.
	1148	IF ( use_syn_turb_gen ) THEN
[3182]	1149	IF ( bc_dirichlet_l ) THEN
[2938]	1150	id_stg_left_l = myidx
	1151	ELSE
	1152	id_stg_left_l = 0
	1153	ENDIF
[3182]	1154	IF ( bc_dirichlet_r ) THEN
[2938]	1155	id_stg_right_l = myidx
	1156	ELSE
	1157	id_stg_right_l = 0
	1158	ENDIF
[3182]	1159	IF ( bc_dirichlet_s ) THEN
[2938]	1160	id_stg_south_l = myidy
	1161	ELSE
	1162	id_stg_south_l = 0
	1163	ENDIF
[3182]	1164	IF ( bc_dirichlet_n ) THEN
[2938]	1165	id_stg_north_l = myidy
	1166	ELSE
	1167	id_stg_north_l = 0
	1168	ENDIF
[1968]	1169
[2938]	1170	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
	1171	CALL MPI_ALLREDUCE( id_stg_left_l, id_stg_left, 1, MPI_INTEGER, &
	1172	MPI_SUM, comm1dx, ierr )
	1173
	1174	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
	1175	CALL MPI_ALLREDUCE( id_stg_right_l, id_stg_right, 1, MPI_INTEGER, &
	1176	MPI_SUM, comm1dx, ierr )
	1177
	1178	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
	1179	CALL MPI_ALLREDUCE( id_stg_south_l, id_stg_south, 1, MPI_INTEGER, &
	1180	MPI_SUM, comm1dy, ierr )
	1181
	1182	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
	1183	CALL MPI_ALLREDUCE( id_stg_north_l, id_stg_north, 1, MPI_INTEGER, &
	1184	MPI_SUM, comm1dy, ierr )
	1185
	1186	ENDIF
	1187
[151]	1188	!
	1189	!-- Broadcast the id of the inflow PE
[3182]	1190	IF ( bc_dirichlet_l ) THEN
[163]	1191	id_inflow_l = myidx
[151]	1192	ELSE
	1193	id_inflow_l = 0
	1194	ENDIF
[622]	1195	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[151]	1196	CALL MPI_ALLREDUCE( id_inflow_l, id_inflow, 1, MPI_INTEGER, MPI_SUM, &
	1197	comm1dx, ierr )
	1198
[163]	1199	!
	1200	!-- Broadcast the id of the recycling plane
	1201	!-- WARNING: needs to be adjusted in case of inflows other than from left side!
[3884]	1202	IF ( turbulent_inflow ) THEN
	1203
[4045]	1204	IF ( NINT( recycling_width / dx, KIND=idp ) >= nxl .AND. &
	1205	NINT( recycling_width / dx, KIND=idp ) <= nxr ) THEN
[3884]	1206	id_recycling_l = myidx
	1207	ELSE
	1208	id_recycling_l = 0
	1209	ENDIF
	1210	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
	1211	CALL MPI_ALLREDUCE( id_recycling_l, id_recycling, 1, MPI_INTEGER, MPI_SUM, &
	1212	comm1dx, ierr )
	1213
[163]	1214	ENDIF
	1215
[2050]	1216	!
	1217	!-- Broadcast the id of the outflow PE and outflow-source plane
	1218	IF ( turbulent_outflow ) THEN
	1219
[3182]	1220	IF ( bc_radiation_r ) THEN
[2050]	1221	id_outflow_l = myidx
	1222	ELSE
	1223	id_outflow_l = 0
	1224	ENDIF
	1225	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
	1226	CALL MPI_ALLREDUCE( id_outflow_l, id_outflow, 1, MPI_INTEGER, MPI_SUM, &
	1227	comm1dx, ierr )
	1228
	1229	IF ( NINT( outflow_source_plane / dx ) >= nxl .AND. &
	1230	NINT( outflow_source_plane / dx ) <= nxr ) THEN
	1231	id_outflow_source_l = myidx
	1232	ELSE
	1233	id_outflow_source_l = 0
	1234	ENDIF
	1235	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
	1236	CALL MPI_ALLREDUCE( id_outflow_source_l, id_outflow_source, 1, &
	1237	MPI_INTEGER, MPI_SUM, comm1dx, ierr )
	1238
	1239	ENDIF
	1240
[3885]	1241	CALL location_message( 'creating virtual PE grids + MPI derived data types', 'finished' )
[1384]	1242
[1804]	1243	#else
[1159]	1244	IF ( bc_lr == 'dirichlet/radiation' ) THEN
[3182]	1245	bc_dirichlet_l = .TRUE.
	1246	bc_radiation_r = .TRUE.
[1159]	1247	ELSEIF ( bc_lr == 'radiation/dirichlet' ) THEN
[3182]	1248	bc_radiation_l = .TRUE.
	1249	bc_dirichlet_r = .TRUE.
[1]	1250	ENDIF
	1251
[1159]	1252	IF ( bc_ns == 'dirichlet/radiation' ) THEN
[3182]	1253	bc_dirichlet_n = .TRUE.
	1254	bc_radiation_s = .TRUE.
[1159]	1255	ELSEIF ( bc_ns == 'radiation/dirichlet' ) THEN
[3182]	1256	bc_radiation_n = .TRUE.
	1257	bc_dirichlet_s = .TRUE.
[1]	1258	ENDIF
	1259	#endif
[807]	1260
[106]	1261	!
[978]	1262	!-- At the inflow or outflow, u or v, respectively, have to be calculated for
	1263	!-- one more grid point.
[3182]	1264	IF ( bc_dirichlet_l .OR. bc_radiation_l ) THEN
[106]	1265	nxlu = nxl + 1
	1266	ELSE
	1267	nxlu = nxl
	1268	ENDIF
[3182]	1269	IF ( bc_dirichlet_s .OR. bc_radiation_s ) THEN
[106]	1270	nysv = nys + 1
	1271	ELSE
	1272	nysv = nys
	1273	ENDIF
[1]	1274
[114]	1275	!
	1276	!-- Allocate wall flag arrays used in the multigrid solver
[1575]	1277	IF ( psolver(1:9) == 'multigrid' ) THEN
[114]	1278
	1279	DO i = maximum_grid_level, 1, -1
	1280
	1281	SELECT CASE ( i )
	1282
	1283	CASE ( 1 )
	1284	ALLOCATE( wall_flags_1(nzb:nzt_mg(i)+1, &
	1285	nys_mg(i)-1:nyn_mg(i)+1, &
	1286	nxl_mg(i)-1:nxr_mg(i)+1) )
	1287
	1288	CASE ( 2 )
	1289	ALLOCATE( wall_flags_2(nzb:nzt_mg(i)+1, &
	1290	nys_mg(i)-1:nyn_mg(i)+1, &
	1291	nxl_mg(i)-1:nxr_mg(i)+1) )
	1292
	1293	CASE ( 3 )
	1294	ALLOCATE( wall_flags_3(nzb:nzt_mg(i)+1, &
	1295	nys_mg(i)-1:nyn_mg(i)+1, &
	1296	nxl_mg(i)-1:nxr_mg(i)+1) )
	1297
	1298	CASE ( 4 )
	1299	ALLOCATE( wall_flags_4(nzb:nzt_mg(i)+1, &
	1300	nys_mg(i)-1:nyn_mg(i)+1, &
	1301	nxl_mg(i)-1:nxr_mg(i)+1) )
	1302
	1303	CASE ( 5 )
	1304	ALLOCATE( wall_flags_5(nzb:nzt_mg(i)+1, &
	1305	nys_mg(i)-1:nyn_mg(i)+1, &
	1306	nxl_mg(i)-1:nxr_mg(i)+1) )
	1307
	1308	CASE ( 6 )
	1309	ALLOCATE( wall_flags_6(nzb:nzt_mg(i)+1, &
	1310	nys_mg(i)-1:nyn_mg(i)+1, &
	1311	nxl_mg(i)-1:nxr_mg(i)+1) )
	1312
	1313	CASE ( 7 )
	1314	ALLOCATE( wall_flags_7(nzb:nzt_mg(i)+1, &
	1315	nys_mg(i)-1:nyn_mg(i)+1, &
	1316	nxl_mg(i)-1:nxr_mg(i)+1) )
	1317
	1318	CASE ( 8 )
	1319	ALLOCATE( wall_flags_8(nzb:nzt_mg(i)+1, &
	1320	nys_mg(i)-1:nyn_mg(i)+1, &
	1321	nxl_mg(i)-1:nxr_mg(i)+1) )
	1322
	1323	CASE ( 9 )
	1324	ALLOCATE( wall_flags_9(nzb:nzt_mg(i)+1, &
	1325	nys_mg(i)-1:nyn_mg(i)+1, &
	1326	nxl_mg(i)-1:nxr_mg(i)+1) )
	1327
	1328	CASE ( 10 )
	1329	ALLOCATE( wall_flags_10(nzb:nzt_mg(i)+1, &
	1330	nys_mg(i)-1:nyn_mg(i)+1, &
	1331	nxl_mg(i)-1:nxr_mg(i)+1) )
	1332
	1333	CASE DEFAULT
[254]	1334	message_string = 'more than 10 multigrid levels'
	1335	CALL message( 'init_pegrid', 'PA0238', 1, 2, 0, 6, 0 )
[114]	1336
	1337	END SELECT
	1338
	1339	ENDDO
	1340
	1341	ENDIF
	1342
[1]	1343	END SUBROUTINE init_pegrid

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |