Home

Context Navigation

poismg.f90 @ 86

Last change on this file since 86 was 77, checked in by raasch, 18 years ago

New:
---

particle reflection from vertical walls implemented, particle SGS model adjusted to walls

Wall functions for vertical walls now include diabatic conditions. New subroutines wall_fluxes, wall_fluxes_e. New 4D-array rif_wall.

new d3par-parameter netcdf_64bit_3d to switch on 64bit offset only for 3D files

new d3par-parameter dt_max to define the maximum value for the allowed timestep

new inipar-parameter loop_optimization to control the loop optimization method

new inipar-parameter pt_refrence. If given, this value is used as the reference that in buoyancy terms (otherwise, the instantaneous horizontally averaged temperature is used).

new user interface user_advec_particles

new initializing action "by_user" calls user_init_3d_model and allows the initial setting of all 3d arrays

topography height informations are stored on arrays zu_s_inner and zw_w_inner and output to the 2d/3d NetCDF files

samples added to the user interface which show how to add user-define time series quantities.

calculation/output of precipitation amount, precipitation rate and z0 (by setting "pra*", "prr*", "z0*" with data_output). The time interval on which the precipitation amount is defined is set by new d3par-parameter precipitation_amount_interval

unit 9 opened for debug output (file DEBUG_<pe#>)

Makefile, advec_particles, average_3d_data, buoyancy, calc_precipitation, check_open, check_parameters, data_output_2d, diffusion_e, diffusion_u, diffusion_v, diffusion_w, diffusivities, header, impact_of_latent_heat, init_particles, init_3d_model, modules, netcdf, parin, production_e, read_var_list, read_3d_binary, sum_up_3d_data, user_interface, write_var_list, write_3d_binary

New: wall_fluxes

Changed:

General revision of non-cyclic horizontal boundary conditions: radiation boundary conditions are now used instead of Neumann conditions at the outflow (calculation needs velocity values for t-dt, which are stored on new arrays u_m_l, u_m_r, etc.), calculation of mean outflow is not needed any more, volume flow control is added for the outflow boundary (currently only for the north boundary!!), additional gridpoints along x and y (uxrp, vynp) are not needed any more, routine "boundary_conds" now operates on timelevel t+dt and is not split in two parts (main, uvw_outflow) any more, Neumann boundary conditions at inflow/outflow in case of non-cyclic boundary conditions for all 2d-arrays that are handled by exchange_horiz_2d

The FFT-method for solving the Poisson-equation is now working with Neumann boundary conditions both at the bottom and the top. This requires adjustments of the tridiagonal coefficients and subtracting the horizontally averaged mean from the vertical velocity field.

+age_m in particle_type

Particles-package is now part of the default code ("-p particles" is not needed any more).

Move call of user_actions( 'after_integration' ) below increment of times
and counters. user_actions is now called for each statistic region and has as an argument the number of the respective region (sr)

d3par-parameter data_output_ts removed. Timeseries output for "profil" removed. Timeseries are now switched on by dt_dots. Timeseries data is collected in flow_statistics.

Initial velocities at nzb+1 are regarded for volume flow control in case they have been set zero before (to avoid small timesteps); see new internal parameters u/v_nzb_p1_for_vfc.

q is not allowed to become negative (prognostic_equations).

poisfft_init is only called if fft-solver is switched on (init_pegrid).

d3par-parameter moisture renamed to humidity.

Subversion global revision number is read from mrun and added to the run description header and to the run control (_rc) file.

vtk directives removed from main program.

The uitility routine interpret_config reads PALM environment variables from NAMELIST instead using the system call GETENV.

advec_u_pw, advec_u_up, advec_v_pw, advec_v_up, asselin_filter, check_parameters, coriolis, data_output_dvrp, data_output_ptseries, data_output_ts, data_output_2d, data_output_3d, diffusion_u, diffusion_v, exchange_horiz, exchange_horiz_2d, flow_statistics, header, init_grid, init_particles, init_pegrid, init_rankine, init_pt_anomaly, init_1d_model, init_3d_model, modules, palm, package_parin, parin, poisfft, poismg, prandtl_fluxes, pres, production_e, prognostic_equations, read_var_list, read_3d_binary, sor, swap_timelevel, time_integration, write_var_list, write_3d_binary

Errors:

Bugfix: preset of tendencies te_em, te_um, te_vm in init_1d_model

Bugfix in sample for reading user defined data from restart file (user_init)

Bugfix in setting diffusivities for cases with the outflow damping layer extending over more than one subdomain (init_3d_model)

Check for possible negative humidities in the initial humidity profile.

in Makefile, default suffixes removed from the suffix list to avoid calling of m2c in
# case of .mod files

Makefile
check_parameters, init_1d_model, init_3d_model, user_interface

Property svn:keywords set to Id

File size: 39.5 KB

Rev	Line
[1]	1	SUBROUTINE poismg( r )
	2
	3	!------------------------------------------------------------------------------!
	4	! Attention: Loop unrolling and cache optimization in SOR-Red/Black method
	5	! still does not bring the expected speedup on ibm! Further work
	6	! is required.
	7	!
	8	! Actual revisions:
	9	! -----------------
[77]	10	!
[1]	11	!
	12	! Former revisions:
	13	! -----------------
[3]	14	! $Id: poismg.f90 77 2007-03-29 04:26:56Z raasch $
[77]	15	!
	16	! 75 2007-03-22 09:54:05Z raasch
	17	! 2nd+3rd argument removed from exchange horiz
	18	!
[3]	19	! RCS Log replace by Id keyword, revision history cleaned up
	20	!
[1]	21	! Revision 1.6 2005/03/26 20:55:54 raasch
	22	! Implementation of non-cyclic (Neumann) horizontal boundary conditions,
	23	! routine prolong simplified (one call of exchange_horiz spared)
	24	!
	25	! Revision 1.1 2001/07/20 13:10:51 raasch
	26	! Initial revision
	27	!
	28	!
	29	! Description:
	30	! ------------
	31	! Solves the Poisson equation for the perturbation pressure with a multigrid
	32	! V- or W-Cycle scheme.
	33	!
	34	! This multigrid method was originally developed for PALM by Joerg Uhlenbrock,
	35	! September 2000 - July 2001.
	36	!------------------------------------------------------------------------------!
	37
	38	USE arrays_3d
	39	USE control_parameters
	40	USE cpulog
	41	USE grid_variables
	42	USE indices
	43	USE interfaces
	44	USE pegrid
	45
	46	IMPLICIT NONE
	47
	48	REAL :: maxerror, maximum_mgcycles, residual_norm
	49
	50	REAL, DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) :: r
	51
	52	REAL, DIMENSION(:,:,:), ALLOCATABLE :: p3
	53
	54
	55	CALL cpu_log( log_point_s(29), 'poismg', 'start' )
	56
	57
	58	!
	59	!-- Initialize arrays and variables used in this subroutine
	60	ALLOCATE ( p3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
	61
	62
	63	!
	64	!-- Some boundaries have to be added to divergence array
[75]	65	CALL exchange_horiz( d )
[1]	66	d(nzb,:,:) = d(nzb+1,:,:)
	67
	68	!
	69	!-- Initiation of the multigrid scheme. Does n cycles until the
	70	!-- residual is smaller than the given limit. The accuracy of the solution
	71	!-- of the poisson equation will increase with the number of cycles.
	72	!-- If the number of cycles is preset by the user, this number will be
	73	!-- carried out regardless of the accuracy.
	74	grid_level_count = 0
	75	mgcycles = 0
	76	IF ( mg_cycles == -1 ) THEN
	77	maximum_mgcycles = 0
	78	residual_norm = 1.0
	79	ELSE
	80	maximum_mgcycles = mg_cycles
	81	residual_norm = 0.0
	82	ENDIF
	83
	84	DO WHILE ( residual_norm > residual_limit .OR. &
	85	mgcycles < maximum_mgcycles )
	86
	87	CALL next_mg_level( d, p, p3, r)
	88
	89	!
	90	!-- Calculate the residual if the user has not preset the number of
	91	!-- cycles to be performed
	92	IF ( maximum_mgcycles == 0 ) THEN
	93	CALL resid( d, p, r )
	94	maxerror = SUM( r(nzb+1:nzt,nys:nyn,nxl:nxr)**2 )
	95	#if defined( __parallel )
	96	CALL MPI_ALLREDUCE( maxerror, residual_norm, 1, MPI_REAL, MPI_SUM, &
	97	comm2d, ierr)
	98	#else
	99	residual_norm = maxerror
	100	#endif
	101	residual_norm = SQRT( residual_norm )
	102	ENDIF
	103
	104	mgcycles = mgcycles + 1
	105
	106	!
	107	!-- If the user has not limited the number of cycles, stop the run in case
	108	!-- of insufficient convergence
	109	IF ( mgcycles > 1000 .AND. mg_cycles == -1 ) THEN
	110	IF ( myid == 0 ) THEN
	111	PRINT*, '+++ poismg: no sufficient convergence within 1000 cycles'
	112	ENDIF
	113	CALL local_stop
	114	ENDIF
	115
	116	ENDDO
	117
	118	DEALLOCATE( p3 )
	119
	120	CALL cpu_log( log_point_s(29), 'poismg', 'stop' )
	121
	122	END SUBROUTINE poismg
	123
	124
	125
	126	SUBROUTINE resid( f_mg, p_mg, r )
	127
	128	!------------------------------------------------------------------------------!
	129	! Description:
	130	! ------------
	131	! Computes the residual of the perturbation pressure.
	132	!------------------------------------------------------------------------------!
	133
	134	USE arrays_3d
	135	USE control_parameters
	136	USE grid_variables
	137	USE indices
	138	USE pegrid
	139
	140	IMPLICIT NONE
	141
	142	INTEGER :: i, j, k, l
	143
	144	REAL, DIMENSION(nzb:nzt_mg(grid_level)+1, &
	145	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
	146	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg, p_mg, r
	147
	148	!
	149	!-- Calculate the residual
	150	l = grid_level
	151
	152	!$OMP PARALLEL PRIVATE (i,j,k)
	153	!$OMP DO
	154	DO i = nxl_mg(l), nxr_mg(l)
	155	DO j = nys_mg(l), nyn_mg(l)
	156	DO k = nzb+1, nzt_mg(l)
	157	r(k,j,i) = f_mg(k,j,i) &
	158	- ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
	159	- ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
	160	- f2_mg(k,l) * p_mg(k+1,j,i) &
	161	- f3_mg(k,l) * p_mg(k-1,j,i) &
	162	+ f1_mg(k,l) * p_mg(k,j,i)
	163	ENDDO
	164	ENDDO
	165	ENDDO
	166	!$OMP END PARALLEL
	167
	168	!
	169	!-- Horizontal boundary conditions
[75]	170	CALL exchange_horiz( r )
[1]	171
	172	IF ( bc_lr /= 'cyclic' ) THEN
	173	IF ( inflow_l .OR. outflow_l ) r(:,:,nxl_mg(l)-1) = r(:,:,nxl_mg(l))
	174	IF ( inflow_r .OR. outflow_r ) r(:,:,nxr_mg(l)+1) = r(:,:,nxr_mg(l))
	175	ENDIF
	176
	177	IF ( bc_ns /= 'cyclic' ) THEN
	178	IF ( inflow_n .OR. outflow_n ) r(:,nyn_mg(l)+1,:) = r(:,nyn_mg(l),:)
	179	IF ( inflow_s .OR. outflow_s ) r(:,nys_mg(l)-1,:) = r(:,nys_mg(l),:)
	180	ENDIF
	181
	182	!
	183	!-- Bottom and top boundary conditions
	184	IF ( ibc_p_b == 1 ) THEN
	185	r(nzb,:,: ) = r(nzb+1,:,:)
	186	ELSE
	187	r(nzb,:,: ) = 0.0
	188	ENDIF
	189
	190	IF ( ibc_p_t == 1 ) THEN
	191	r(nzt_mg(l)+1,:,: ) = r(nzt_mg(l),:,:)
	192	ELSE
	193	r(nzt_mg(l)+1,:,: ) = 0.0
	194	ENDIF
	195
	196
	197	END SUBROUTINE resid
	198
	199
	200
	201	SUBROUTINE restrict( f_mg, r )
	202
	203	!------------------------------------------------------------------------------!
	204	! Description:
	205	! ------------
	206	! Interpolates the residual on the next coarser grid with "full weighting"
	207	! scheme
	208	!------------------------------------------------------------------------------!
	209
	210	USE control_parameters
	211	USE grid_variables
	212	USE indices
	213	USE pegrid
	214
	215	IMPLICIT NONE
	216
	217	INTEGER :: i, ic, j, jc, k, kc, l
	218
	219	REAL, DIMENSION(nzb:nzt_mg(grid_level)+1, &
	220	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
	221	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg
	222
	223	REAL, DIMENSION(nzb:nzt_mg(grid_level+1)+1, &
	224	nys_mg(grid_level+1)-1:nyn_mg(grid_level+1)+1, &
	225	nxl_mg(grid_level+1)-1:nxr_mg(grid_level+1)+1) :: r
	226
	227	!
	228	!-- Interpolate the residual
	229	l = grid_level
	230
	231	!$OMP PARALLEL PRIVATE (i,j,k,ic,jc,kc)
	232	!$OMP DO
	233	DO ic = nxl_mg(l), nxr_mg(l)
	234	i = 2*ic
	235	DO jc = nys_mg(l), nyn_mg(l)
	236	j = 2*jc
	237	DO kc = nzb+1, nzt_mg(l)
	238	k = 2*kc-1
	239	f_mg(kc,jc,ic) = 1.0 / 64.0 * ( &
	240	8.0 * r(k,j,i) &
	241	+ 4.0 * ( r(k,j,i-1) + r(k,j,i+1) + &
	242	r(k,j+1,i) + r(k,j-1,i) ) &
	243	+ 2.0 * ( r(k,j-1,i-1) + r(k,j+1,i-1) + &
	244	r(k,j-1,i+1) + r(k,j+1,i+1) ) &
	245	+ 4.0 * r(k-1,j,i) &
	246	+ 2.0 * ( r(k-1,j,i-1) + r(k-1,j,i+1) + &
	247	r(k-1,j+1,i) + r(k-1,j-1,i) ) &
	248	+ ( r(k-1,j-1,i-1) + r(k-1,j+1,i-1) + &
	249	r(k-1,j-1,i+1) + r(k-1,j+1,i+1) ) &
	250	+ 4.0 * r(k+1,j,i) &
	251	+ 2.0 * ( r(k+1,j,i-1) + r(k+1,j,i+1) + &
	252	r(k+1,j+1,i) + r(k+1,j-1,i) ) &
	253	+ ( r(k+1,j-1,i-1) + r(k+1,j+1,i-1) + &
	254	r(k+1,j-1,i+1) + r(k+1,j+1,i+1) ) &
	255	)
	256	ENDDO
	257	ENDDO
	258	ENDDO
	259	!$OMP END PARALLEL
	260
	261	!
	262	!-- Horizontal boundary conditions
[75]	263	CALL exchange_horiz( f_mg )
[1]	264
	265	IF ( bc_lr /= 'cyclic' ) THEN
	266	IF (inflow_l .OR. outflow_l) f_mg(:,:,nxl_mg(l)-1) = f_mg(:,:,nxl_mg(l))
	267	IF (inflow_r .OR. outflow_r) f_mg(:,:,nxr_mg(l)+1) = f_mg(:,:,nxr_mg(l))
	268	ENDIF
	269
	270	IF ( bc_ns /= 'cyclic' ) THEN
	271	IF (inflow_n .OR. outflow_n) f_mg(:,nyn_mg(l)+1,:) = f_mg(:,nyn_mg(l),:)
	272	IF (inflow_s .OR. outflow_s) f_mg(:,nys_mg(l)-1,:) = f_mg(:,nys_mg(l),:)
	273	ENDIF
	274
	275	!
	276	!-- Bottom and top boundary conditions
	277	IF ( ibc_p_b == 1 ) THEN
	278	f_mg(nzb,:,: ) = f_mg(nzb+1,:,:)
	279	ELSE
	280	f_mg(nzb,:,: ) = 0.0
	281	ENDIF
	282
	283	IF ( ibc_p_t == 1 ) THEN
	284	f_mg(nzt_mg(l)+1,:,: ) = f_mg(nzt_mg(l),:,:)
	285	ELSE
	286	f_mg(nzt_mg(l)+1,:,: ) = 0.0
	287	ENDIF
	288
	289
	290	END SUBROUTINE restrict
	291
	292
	293
	294	SUBROUTINE prolong( p, temp )
	295
	296	!------------------------------------------------------------------------------!
	297	! Description:
	298	! ------------
	299	! Interpolates the correction of the perturbation pressure
	300	! to the next finer grid.
	301	!------------------------------------------------------------------------------!
	302
	303	USE control_parameters
	304	USE pegrid
	305	USE indices
	306
	307	IMPLICIT NONE
	308
	309	INTEGER :: i, j, k, l
	310
	311	REAL, DIMENSION(nzb:nzt_mg(grid_level-1)+1, &
	312	nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1, &
	313	nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1 ) :: p
	314
	315	REAL, DIMENSION(nzb:nzt_mg(grid_level)+1, &
	316	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
	317	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: temp
	318
	319
	320	!
	321	!-- First, store elements of the coarser grid on the next finer grid
	322	l = grid_level
	323
	324	!$OMP PARALLEL PRIVATE (i,j,k)
	325	!$OMP DO
	326	DO i = nxl_mg(l-1), nxr_mg(l-1)
	327	DO j = nys_mg(l-1), nyn_mg(l-1)
	328	!CDIR NODEP
	329	DO k = nzb+1, nzt_mg(l-1)
	330	!
	331	!-- Points of the coarse grid are directly stored on the next finer
	332	!-- grid
	333	temp(2k-1,2j,2*i) = p(k,j,i)
	334	!
	335	!-- Points between two coarse-grid points
	336	temp(2k-1,2j,2i+1) = 0.5 ( p(k,j,i) + p(k,j,i+1) )
	337	temp(2k-1,2j+1,2i) = 0.5 ( p(k,j,i) + p(k,j+1,i) )
	338	temp(2k,2j,2i) = 0.5 ( p(k,j,i) + p(k+1,j,i) )
	339	!
	340	!-- Points in the center of the planes stretched by four points
	341	!-- of the coarse grid cube
	342	temp(2k-1,2j+1,2i+1) = 0.25 ( p(k,j,i) + p(k,j,i+1) + &
	343	p(k,j+1,i) + p(k,j+1,i+1) )
	344	temp(2k,2j,2i+1) = 0.25 ( p(k,j,i) + p(k,j,i+1) + &
	345	p(k+1,j,i) + p(k+1,j,i+1) )
	346	temp(2k,2j+1,2i) = 0.25 ( p(k,j,i) + p(k,j+1,i) + &
	347	p(k+1,j,i) + p(k+1,j+1,i) )
	348	!
	349	!-- Points in the middle of coarse grid cube
	350	temp(2k,2j+1,2i+1) = 0.125 ( p(k,j,i) + p(k,j,i+1) + &
	351	p(k,j+1,i) + p(k,j+1,i+1) + &
	352	p(k+1,j,i) + p(k+1,j,i+1) + &
	353	p(k+1,j+1,i) + p(k+1,j+1,i+1) )
	354	ENDDO
	355	ENDDO
	356	ENDDO
	357	!$OMP END PARALLEL
	358
	359	!
	360	!-- Horizontal boundary conditions
[75]	361	CALL exchange_horiz( temp )
[1]	362
	363	IF ( bc_lr /= 'cyclic' ) THEN
	364	IF (inflow_l .OR. outflow_l) temp(:,:,nxl_mg(l)-1) = temp(:,:,nxl_mg(l))
	365	IF (inflow_r .OR. outflow_r) temp(:,:,nxr_mg(l)+1) = temp(:,:,nxr_mg(l))
	366	ENDIF
	367
	368	IF ( bc_ns /= 'cyclic' ) THEN
	369	IF (inflow_n .OR. outflow_n) temp(:,nyn_mg(l)+1,:) = temp(:,nyn_mg(l),:)
	370	IF (inflow_s .OR. outflow_s) temp(:,nys_mg(l)-1,:) = temp(:,nys_mg(l),:)
	371	ENDIF
	372
	373	!
	374	!-- Bottom and top boundary conditions
	375	IF ( ibc_p_b == 1 ) THEN
	376	temp(nzb,:,: ) = temp(nzb+1,:,:)
	377	ELSE
	378	temp(nzb,:,: ) = 0.0
	379	ENDIF
	380
	381	IF ( ibc_p_t == 1 ) THEN
	382	temp(nzt_mg(l)+1,:,: ) = temp(nzt_mg(l),:,:)
	383	ELSE
	384	temp(nzt_mg(l)+1,:,: ) = 0.0
	385	ENDIF
	386
	387
	388	END SUBROUTINE prolong
	389
	390
	391	SUBROUTINE redblack( f_mg, p_mg )
	392
	393	!------------------------------------------------------------------------------!
	394	! Description:
	395	! ------------
	396	! Relaxation method for the multigrid scheme. A Gauss-Seidel iteration with
	397	! 3D-Red-Black decomposition (GS-RB) is used.
	398	!------------------------------------------------------------------------------!
	399
	400	USE arrays_3d
	401	USE control_parameters
	402	USE cpulog
	403	USE grid_variables
	404	USE indices
	405	USE interfaces
	406	USE pegrid
	407
	408	IMPLICIT NONE
	409
	410	INTEGER :: colour, i, ic, j, jc, jj, k, l, n
	411
	412	LOGICAL :: unroll
	413
	414	REAL, DIMENSION(nzb:nzt_mg(grid_level)+1, &
	415	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
	416	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg, p_mg
	417
	418
	419	l = grid_level
	420
	421	unroll = ( MOD( nyn_mg(l)-nys_mg(l)+1, 4 ) == 0 .AND. &
	422	MOD( nxr_mg(l)-nxl_mg(l)+1, 2 ) == 0 )
	423
	424	DO n = 1, ngsrb
	425
	426	DO colour = 1, 2
	427
	428	IF ( .NOT. unroll ) THEN
	429	CALL cpu_log( log_point_s(36), 'redblack_no_unroll', 'start' )
	430
	431	!
	432	!-- Without unrolling of loops, no cache optimization
	433	DO i = nxl_mg(l), nxr_mg(l), 2
	434	DO j = nys_mg(l) + 2 - colour, nyn_mg(l), 2
	435	DO k = nzb+1, nzt_mg(l), 2
	436	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
	437	ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
	438	+ ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
	439	+ f2_mg(k,l) * p_mg(k+1,j,i) &
	440	+ f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
	441	)
	442	ENDDO
	443	ENDDO
	444	ENDDO
	445
	446	DO i = nxl_mg(l)+1, nxr_mg(l), 2
	447	DO j = nys_mg(l) + (colour-1), nyn_mg(l), 2
	448	DO k = nzb+1, nzt_mg(l), 2
	449	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
	450	ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
	451	+ ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
	452	+ f2_mg(k,l) * p_mg(k+1,j,i) &
	453	+ f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
	454	)
	455	ENDDO
	456	ENDDO
	457	ENDDO
	458
	459	DO i = nxl_mg(l), nxr_mg(l), 2
	460	DO j = nys_mg(l) + (colour-1), nyn_mg(l), 2
	461	DO k = nzb+2, nzt_mg(l), 2
	462	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
	463	ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
	464	+ ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
	465	+ f2_mg(k,l) * p_mg(k+1,j,i) &
	466	+ f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
	467	)
	468	ENDDO
	469	ENDDO
	470	ENDDO
	471
	472	DO i = nxl_mg(l)+1, nxr_mg(l), 2
	473	DO j = nys_mg(l) + 2 - colour, nyn_mg(l), 2
	474	DO k = nzb+2, nzt_mg(l), 2
	475	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
	476	ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
	477	+ ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
	478	+ f2_mg(k,l) * p_mg(k+1,j,i) &
	479	+ f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
	480	)
	481	ENDDO
	482	ENDDO
	483	ENDDO
	484	CALL cpu_log( log_point_s(36), 'redblack_no_unroll', 'stop' )
	485
	486	ELSE
	487
	488	!
	489	!-- Loop unrolling along y, only one i loop for better cache use
	490	CALL cpu_log( log_point_s(38), 'redblack_unroll', 'start' )
	491	DO ic = nxl_mg(l), nxr_mg(l), 2
	492	DO jc = nys_mg(l), nyn_mg(l), 4
	493	i = ic
	494	jj = jc+2-colour
	495	DO k = nzb+1, nzt_mg(l), 2
	496	j = jj
	497	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
	498	ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
	499	+ ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
	500	+ f2_mg(k,l) * p_mg(k+1,j,i) &
	501	+ f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
	502	)
	503	j = jj+2
	504	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
	505	ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
	506	+ ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
	507	+ f2_mg(k,l) * p_mg(k+1,j,i) &
	508	+ f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
	509	)
	510	! j = jj+4
	511	! p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
	512	! ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
	513	! + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
	514	! + f2_mg(k,l) * p_mg(k+1,j,i) &
	515	! + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
	516	! )
	517	! j = jj+6
	518	! p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
	519	! ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
	520	! + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
	521	! + f2_mg(k,l) * p_mg(k+1,j,i) &
	522	! + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
	523	! )
	524	ENDDO
	525
	526	i = ic+1
	527	jj = jc+colour-1
	528	DO k = nzb+1, nzt_mg(l), 2
	529	j =jj
	530	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
	531	ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
	532	+ ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
	533	+ f2_mg(k,l) * p_mg(k+1,j,i) &
	534	+ f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
	535	)
	536	j = jj+2
	537	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
	538	ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
	539	+ ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
	540	+ f2_mg(k,l) * p_mg(k+1,j,i) &
	541	+ f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
	542	)
	543	! j = jj+4
	544	! p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
	545	! ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
	546	! + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
	547	! + f2_mg(k,l) * p_mg(k+1,j,i) &
	548	! + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
	549	! )
	550	! j = jj+6
	551	! p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
	552	! ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
	553	! + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
	554	! + f2_mg(k,l) * p_mg(k+1,j,i) &
	555	! + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
	556	! )
	557	ENDDO
	558
	559	i = ic
	560	jj = jc+colour-1
	561	DO k = nzb+2, nzt_mg(l), 2
	562	j =jj
	563	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
	564	ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
	565	+ ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
	566	+ f2_mg(k,l) * p_mg(k+1,j,i) &
	567	+ f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
	568	)
	569	j = jj+2
	570	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
	571	ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
	572	+ ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
	573	+ f2_mg(k,l) * p_mg(k+1,j,i) &
	574	+ f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
	575	)
	576	! j = jj+4
	577	! p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
	578	! ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
	579	! + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
	580	! + f2_mg(k,l) * p_mg(k+1,j,i) &
	581	! + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
	582	! )
	583	! j = jj+6
	584	! p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
	585	! ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
	586	! + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
	587	! + f2_mg(k,l) * p_mg(k+1,j,i) &
	588	! + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
	589	! )
	590	ENDDO
	591
	592	i = ic+1
	593	jj = jc+2-colour
	594	DO k = nzb+2, nzt_mg(l), 2
	595	j =jj
	596	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
	597	ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
	598	+ ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
	599	+ f2_mg(k,l) * p_mg(k+1,j,i) &
	600	+ f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
	601	)
	602	j = jj+2
	603	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
	604	ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
	605	+ ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
	606	+ f2_mg(k,l) * p_mg(k+1,j,i) &
	607	+ f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
	608	)
	609	! j = jj+4
	610	! p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
	611	! ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
	612	! + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
	613	! + f2_mg(k,l) * p_mg(k+1,j,i) &
	614	! + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
	615	! )
	616	! j = jj+6
	617	! p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
	618	! ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
	619	! + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
	620	! + f2_mg(k,l) * p_mg(k+1,j,i) &
	621	! + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
	622	! )
	623	ENDDO
	624
	625	ENDDO
	626	ENDDO
	627	CALL cpu_log( log_point_s(38), 'redblack_unroll', 'stop' )
	628
	629	ENDIF
	630
	631	!
	632	!-- Horizontal boundary conditions
[75]	633	CALL exchange_horiz( p_mg )
[1]	634
	635	IF ( bc_lr /= 'cyclic' ) THEN
	636	IF ( inflow_l .OR. outflow_l ) THEN
	637	p_mg(:,:,nxl_mg(l)-1) = p_mg(:,:,nxl_mg(l))
	638	ENDIF
	639	IF ( inflow_r .OR. outflow_r ) THEN
	640	p_mg(:,:,nxr_mg(l)+1) = p_mg(:,:,nxr_mg(l))
	641	ENDIF
	642	ENDIF
	643
	644	IF ( bc_ns /= 'cyclic' ) THEN
	645	IF ( inflow_n .OR. outflow_n ) THEN
	646	p_mg(:,nyn_mg(l)+1,:) = p_mg(:,nyn_mg(l),:)
	647	ENDIF
	648	IF ( inflow_s .OR. outflow_s ) THEN
	649	p_mg(:,nys_mg(l)-1,:) = p_mg(:,nys_mg(l),:)
	650	ENDIF
	651	ENDIF
	652
	653	!
	654	!-- Bottom and top boundary conditions
	655	IF ( ibc_p_b == 1 ) THEN
	656	p_mg(nzb,:,: ) = p_mg(nzb+1,:,:)
	657	ELSE
	658	p_mg(nzb,:,: ) = 0.0
	659	ENDIF
	660
	661	IF ( ibc_p_t == 1 ) THEN
	662	p_mg(nzt_mg(l)+1,:,: ) = p_mg(nzt_mg(l),:,:)
	663	ELSE
	664	p_mg(nzt_mg(l)+1,:,: ) = 0.0
	665	ENDIF
	666
	667	ENDDO
	668
	669	ENDDO
	670
	671
	672	END SUBROUTINE redblack
	673
	674
	675
	676	SUBROUTINE mg_gather( f2, f2_sub )
	677
	678	USE control_parameters
	679	USE cpulog
	680	USE indices
	681	USE interfaces
	682	USE pegrid
	683
	684	IMPLICIT NONE
	685
	686	INTEGER :: n, nwords, sender
	687
	688	REAL, DIMENSION(nzb:nzt_mg(grid_level)+1, &
	689	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
	690	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f2
	691
	692	REAL, DIMENSION(nzb:mg_loc_ind(5,myid)+1, &
	693	mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, &
	694	mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) :: f2_sub
	695
	696	!
	697	!-- Find out the number of array elements of the subdomain array
	698	nwords = SIZE( f2_sub )
	699
	700	#if defined( __parallel )
	701	CALL cpu_log( log_point_s(34), 'mg_gather', 'start' )
	702
	703	IF ( myid == 0 ) THEN
	704	!
	705	!-- Store the local subdomain array on the total array
	706	f2(:,mg_loc_ind(3,0)-1:mg_loc_ind(4,0)+1, &
	707	mg_loc_ind(1,0)-1:mg_loc_ind(2,0)+1) = f2_sub
	708
	709	!
	710	!-- Receive the subdomain arrays from all other PEs and store them on the
	711	!-- total array
	712	DO n = 1, numprocs-1
	713	!
	714	!-- Receive the arrays in arbitrary order from the PEs.
	715	CALL MPI_RECV( f2_sub(nzb,mg_loc_ind(3,0)-1,mg_loc_ind(1,0)-1), &
	716	nwords, MPI_REAL, MPI_ANY_SOURCE, 1, comm2d, status, &
	717	ierr )
	718	sender = status(MPI_SOURCE)
	719	f2(:,mg_loc_ind(3,sender)-1:mg_loc_ind(4,sender)+1, &
	720	mg_loc_ind(1,sender)-1:mg_loc_ind(2,sender)+1) = f2_sub
	721	ENDDO
	722
	723	ELSE
	724	!
	725	!-- Send subdomain array to PE0
	726	CALL MPI_SEND( f2_sub(nzb,mg_loc_ind(3,myid)-1,mg_loc_ind(1,myid)-1), &
	727	nwords, MPI_REAL, 0, 1, comm2d, ierr )
	728	ENDIF
	729
	730	CALL cpu_log( log_point_s(34), 'mg_gather', 'stop' )
	731	#endif
	732
	733	END SUBROUTINE mg_gather
	734
	735
	736
	737	SUBROUTINE mg_scatter( p2, p2_sub )
	738	!
	739	!-- TODO: It may be possible to improve the speed of this routine by using
	740	!-- non-blocking communication
	741
	742	USE control_parameters
	743	USE cpulog
	744	USE indices
	745	USE interfaces
	746	USE pegrid
	747
	748	IMPLICIT NONE
	749
	750	INTEGER :: n, nwords, sender
	751
	752	REAL, DIMENSION(nzb:nzt_mg(grid_level-1)+1, &
	753	nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1, &
	754	nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) :: p2
	755
	756	REAL, DIMENSION(nzb:mg_loc_ind(5,myid)+1, &
	757	mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, &
	758	mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) :: p2_sub
	759
	760	!
	761	!-- Find out the number of array elements of the subdomain array
	762	nwords = SIZE( p2_sub )
	763
	764	#if defined( __parallel )
	765	CALL cpu_log( log_point_s(35), 'mg_scatter', 'start' )
	766
	767	IF ( myid == 0 ) THEN
	768	!
	769	!-- Scatter the subdomain arrays to the other PEs by blocking
	770	!-- communication
	771	DO n = 1, numprocs-1
	772
	773	p2_sub = p2(:,mg_loc_ind(3,n)-1:mg_loc_ind(4,n)+1, &
	774	mg_loc_ind(1,n)-1:mg_loc_ind(2,n)+1)
	775
	776	CALL MPI_SEND( p2_sub(nzb,mg_loc_ind(3,0)-1,mg_loc_ind(1,0)-1), &
	777	nwords, MPI_REAL, n, 1, comm2d, ierr )
	778
	779	ENDDO
	780
	781	!
	782	!-- Store data from the total array to the local subdomain array
	783	p2_sub = p2(:,mg_loc_ind(3,0)-1:mg_loc_ind(4,0)+1, &
	784	mg_loc_ind(1,0)-1:mg_loc_ind(2,0)+1)
	785
	786	ELSE
	787	!
	788	!-- Receive subdomain array from PE0
	789	CALL MPI_RECV( p2_sub(nzb,mg_loc_ind(3,myid)-1,mg_loc_ind(1,myid)-1), &
	790	nwords, MPI_REAL, 0, 1, comm2d, status, ierr )
	791
	792	ENDIF
	793
	794	CALL cpu_log( log_point_s(35), 'mg_scatter', 'stop' )
	795	#endif
	796
	797	END SUBROUTINE mg_scatter
	798
	799
	800
	801	RECURSIVE SUBROUTINE next_mg_level( f_mg, p_mg, p3, r )
	802
	803	!------------------------------------------------------------------------------!
	804	! Description:
	805	! ------------
	806	! This is where the multigrid technique takes place. V- and W- Cycle are
	807	! implemented and steered by the parameter "gamma". Parameter "nue" determines
	808	! the convergence of the multigrid iterative solution. There are nue times
	809	! RB-GS iterations. It should be set to "1" or "2", considering the time effort
	810	! one would like to invest. Last choice shows a very good converging factor,
	811	! but leads to an increase in computing time.
	812	!------------------------------------------------------------------------------!
	813
	814	USE arrays_3d
	815	USE control_parameters
	816	USE grid_variables
	817	USE indices
	818	USE pegrid
	819
	820	IMPLICIT NONE
	821
	822	INTEGER :: i, j, k, nxl_mg_save, nxr_mg_save, nyn_mg_save, nys_mg_save, &
	823	nzt_mg_save
	824
	825	LOGICAL :: restore_boundary_lr_on_pe0, restore_boundary_ns_on_pe0
	826
	827	REAL, DIMENSION(nzb:nzt_mg(grid_level)+1, &
	828	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
	829	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg, p_mg, p3, r
	830
	831	REAL, DIMENSION(:,:,:), ALLOCATABLE :: f2, f2_sub, p2, p2_sub
	832
	833	!
	834	!-- Restriction to the coarsest grid
	835	10 IF ( grid_level == 1 ) THEN
	836
	837	!
	838	!-- Solution on the coarsest grid. Double the number of Gauss-Seidel
	839	!-- iterations in order to get a more accurate solution.
	840	ngsrb = 2 * ngsrb
	841	CALL redblack( f_mg, p_mg )
	842	ngsrb = ngsrb / 2
	843
	844	ELSEIF ( grid_level /= 1 ) THEN
	845
	846	grid_level_count(grid_level) = grid_level_count(grid_level) + 1
	847
	848	!
	849	!-- Solution on the actual grid level
	850	CALL redblack( f_mg, p_mg )
	851
	852	!
	853	!-- Determination of the actual residual
	854	CALL resid( f_mg, p_mg, r )
	855
	856	!
	857	!-- Restriction of the residual (finer grid values!) to the next coarser
	858	!-- grid. Therefore, the grid level has to be decremented now. nxl..nzt have
	859	!-- to be set to the coarse grid values, because these variables are needed
	860	!-- for the exchange of ghost points in routine exchange_horiz
	861	grid_level = grid_level - 1
	862	nxl = nxl_mg(grid_level)
	863	nxr = nxr_mg(grid_level)
	864	nys = nys_mg(grid_level)
	865	nyn = nyn_mg(grid_level)
	866	nzt = nzt_mg(grid_level)
	867
	868	ALLOCATE( f2(nzb:nzt_mg(grid_level)+1, &
	869	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
	870	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1), &
	871	p2(nzb:nzt_mg(grid_level)+1, &
	872	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
	873	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) )
	874
	875	IF ( grid_level == mg_switch_to_pe0_level ) THEN
	876	! print*, 'myid=',myid, ' restrict and switch to PE0. level=', grid_level
	877	!
	878	!-- From this level on, calculations are done on PE0 only.
	879	!-- First, carry out restriction on the subdomain.
	880	!-- Therefore, indices of the level have to be changed to subdomain values
	881	!-- in between (otherwise, the restrict routine would expect
	882	!-- the gathered array)
	883	nxl_mg_save = nxl_mg(grid_level)
	884	nxr_mg_save = nxr_mg(grid_level)
	885	nys_mg_save = nys_mg(grid_level)
	886	nyn_mg_save = nyn_mg(grid_level)
	887	nzt_mg_save = nzt_mg(grid_level)
	888	nxl_mg(grid_level) = mg_loc_ind(1,myid)
	889	nxr_mg(grid_level) = mg_loc_ind(2,myid)
	890	nys_mg(grid_level) = mg_loc_ind(3,myid)
	891	nyn_mg(grid_level) = mg_loc_ind(4,myid)
	892	nzt_mg(grid_level) = mg_loc_ind(5,myid)
	893	nxl = mg_loc_ind(1,myid)
	894	nxr = mg_loc_ind(2,myid)
	895	nys = mg_loc_ind(3,myid)
	896	nyn = mg_loc_ind(4,myid)
	897	nzt = mg_loc_ind(5,myid)
	898
	899	ALLOCATE( f2_sub(nzb:nzt_mg(grid_level)+1, &
	900	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
	901	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) )
	902
	903	CALL restrict( f2_sub, r )
	904
	905	!
	906	!-- Restore the correct indices of this level
	907	nxl_mg(grid_level) = nxl_mg_save
	908	nxr_mg(grid_level) = nxr_mg_save
	909	nys_mg(grid_level) = nys_mg_save
	910	nyn_mg(grid_level) = nyn_mg_save
	911	nzt_mg(grid_level) = nzt_mg_save
	912	nxl = nxl_mg(grid_level)
	913	nxr = nxr_mg(grid_level)
	914	nys = nys_mg(grid_level)
	915	nyn = nyn_mg(grid_level)
	916	nzt = nzt_mg(grid_level)
	917
	918	!
	919	!-- Gather all arrays from the subdomains on PE0
	920	CALL mg_gather( f2, f2_sub )
	921
	922	!
	923	!-- Set switch for routine exchange_horiz, that no ghostpoint exchange
	924	!-- has to be carried out from now on
	925	mg_switch_to_pe0 = .TRUE.
	926
	927	!
	928	!-- In case of non-cyclic lateral boundary conditions, both in- and
	929	!-- outflow conditions have to be used on PE0 after the switch, because
	930	!-- it then contains the total domain. Due to the virtual processor
	931	!-- grid, before the switch, PE0 can have in-/outflow at the left
	932	!-- and south wall only (or on opposite walls in case of a 1d
	933	!-- decomposition).
	934	restore_boundary_lr_on_pe0 = .FALSE.
	935	restore_boundary_ns_on_pe0 = .FALSE.
	936	IF ( myid == 0 ) THEN
	937	IF ( inflow_l .AND. .NOT. outflow_r ) THEN
	938	outflow_r = .TRUE.
	939	restore_boundary_lr_on_pe0 = .TRUE.
	940	ENDIF
	941	IF ( outflow_l .AND. .NOT. inflow_r ) THEN
	942	inflow_r = .TRUE.
	943	restore_boundary_lr_on_pe0 = .TRUE.
	944	ENDIF
	945	IF ( inflow_s .AND. .NOT. outflow_n ) THEN
	946	outflow_n = .TRUE.
	947	restore_boundary_ns_on_pe0 = .TRUE.
	948	ENDIF
	949	IF ( outflow_s .AND. .NOT. inflow_n ) THEN
	950	inflow_n = .TRUE.
	951	restore_boundary_ns_on_pe0 = .TRUE.
	952	ENDIF
	953	ENDIF
	954
	955	DEALLOCATE( f2_sub )
	956
	957	ELSE
	958
	959	CALL restrict( f2, r )
	960
	961	ENDIF
	962	p2 = 0.0
	963
	964	!
	965	!-- Repeat the same procedure till the coarsest grid is reached
	966	IF ( myid == 0 .OR. grid_level > mg_switch_to_pe0_level ) THEN
	967	CALL next_mg_level( f2, p2, p3, r )
	968	ENDIF
	969
	970	ENDIF
	971
	972	!
	973	!-- Now follows the prolongation
	974	IF ( grid_level >= 2 ) THEN
	975
	976	!
	977	!-- Grid level has to be incremented on the PEs where next_mg_level
	978	!-- has not been called before (normally it is incremented at the end
	979	!-- of next_mg_level)
	980	IF ( myid /= 0 .AND. grid_level == mg_switch_to_pe0_level ) THEN
	981	grid_level = grid_level + 1
	982	nxl = nxl_mg(grid_level)
	983	nxr = nxr_mg(grid_level)
	984	nys = nys_mg(grid_level)
	985	nyn = nyn_mg(grid_level)
	986	nzt = nzt_mg(grid_level)
	987	ENDIF
	988
	989	!
	990	!-- Prolongation of the new residual. The values are transferred
	991	!-- from the coarse to the next finer grid.
	992	IF ( grid_level == mg_switch_to_pe0_level+1 ) THEN
	993	!
	994	!-- At this level, the new residual first has to be scattered from
	995	!-- PE0 to the other PEs
	996	ALLOCATE( p2_sub(nzb:mg_loc_ind(5,myid)+1, &
	997	mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, &
	998	mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) )
	999
	1000	CALL mg_scatter( p2, p2_sub )
	1001
	1002	!
	1003	!-- Therefore, indices of the previous level have to be changed to
	1004	!-- subdomain values in between (otherwise, the prolong routine would
	1005	!-- expect the gathered array)
	1006	nxl_mg_save = nxl_mg(grid_level-1)
	1007	nxr_mg_save = nxr_mg(grid_level-1)
	1008	nys_mg_save = nys_mg(grid_level-1)
	1009	nyn_mg_save = nyn_mg(grid_level-1)
	1010	nzt_mg_save = nzt_mg(grid_level-1)
	1011	nxl_mg(grid_level-1) = mg_loc_ind(1,myid)
	1012	nxr_mg(grid_level-1) = mg_loc_ind(2,myid)
	1013	nys_mg(grid_level-1) = mg_loc_ind(3,myid)
	1014	nyn_mg(grid_level-1) = mg_loc_ind(4,myid)
	1015	nzt_mg(grid_level-1) = mg_loc_ind(5,myid)
	1016
	1017	!
	1018	!-- Set switch for routine exchange_horiz, that ghostpoint exchange
	1019	!-- has to be carried again out from now on
	1020	mg_switch_to_pe0 = .FALSE.
	1021
	1022	!
	1023	!-- In case of non-cyclic lateral boundary conditions, restore the
	1024	!-- in-/outflow conditions on PE0
	1025	IF ( myid == 0 ) THEN
	1026	IF ( restore_boundary_lr_on_pe0 ) THEN
	1027	IF ( inflow_l ) outflow_r = .FALSE.
	1028	IF ( outflow_l ) inflow_r = .FALSE.
	1029	ENDIF
	1030	IF ( restore_boundary_ns_on_pe0 ) THEN
	1031	IF ( inflow_s ) outflow_n = .FALSE.
	1032	IF ( outflow_s ) inflow_n = .FALSE.
	1033	ENDIF
	1034	ENDIF
	1035
	1036	CALL prolong( p2_sub, p3 )
	1037
	1038	!
	1039	!-- Restore the correct indices of the previous level
	1040	nxl_mg(grid_level-1) = nxl_mg_save
	1041	nxr_mg(grid_level-1) = nxr_mg_save
	1042	nys_mg(grid_level-1) = nys_mg_save
	1043	nyn_mg(grid_level-1) = nyn_mg_save
	1044	nzt_mg(grid_level-1) = nzt_mg_save
	1045
	1046	DEALLOCATE( p2_sub )
	1047
	1048	ELSE
	1049
	1050	CALL prolong( p2, p3 )
	1051
	1052	ENDIF
	1053
	1054	!
	1055	!-- Temporary arrays for the actual grid are not needed any more
	1056	DEALLOCATE( p2, f2 )
	1057
	1058	!
	1059	!-- Computation of the new pressure correction. Therefore,
	1060	!-- values from prior grids are added up automatically stage by stage.
	1061	DO i = nxl_mg(grid_level)-1, nxr_mg(grid_level)+1
	1062	DO j = nys_mg(grid_level)-1, nyn_mg(grid_level)+1
	1063	DO k = nzb, nzt_mg(grid_level)+1
	1064	p_mg(k,j,i) = p_mg(k,j,i) + p3(k,j,i)
	1065	ENDDO
	1066	ENDDO
	1067	ENDDO
	1068
	1069	!
	1070	!-- Relaxation of the new solution
	1071	CALL redblack( f_mg, p_mg )
	1072
	1073	ENDIF
	1074
	1075	!
	1076	!-- The following few lines serve the steering of the multigrid scheme
	1077	IF ( grid_level == maximum_grid_level ) THEN
	1078
	1079	GOTO 20
	1080
	1081	ELSEIF ( grid_level /= maximum_grid_level .AND. grid_level /= 1 .AND. &
	1082	grid_level_count(grid_level) /= gamma_mg ) THEN
	1083
	1084	GOTO 10
	1085
	1086	ENDIF
	1087
	1088	!
	1089	!-- Reset counter for the next call of poismg
	1090	grid_level_count(grid_level) = 0
	1091
	1092	!
	1093	!-- Continue with the next finer level. nxl..nzt have to be
	1094	!-- set to the finer grid values, because these variables are needed for the
	1095	!-- exchange of ghost points in routine exchange_horiz
	1096	grid_level = grid_level + 1
	1097	nxl = nxl_mg(grid_level)
	1098	nxr = nxr_mg(grid_level)
	1099	nys = nys_mg(grid_level)
	1100	nyn = nyn_mg(grid_level)
	1101	nzt = nzt_mg(grid_level)
	1102
	1103	20 CONTINUE
	1104
	1105	END SUBROUTINE next_mg_level

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

Context Navigation

source: palm/trunk/SOURCE/poismg.f90 @ 86

Download in other formats: