Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |