Home

Context Navigation

pres.f90 @ 4899

Last change on this file since 4899 was 4855, checked in by raasch, 4 years ago
bugfix: mean w removal not applied to ghost points of the total domain in case of non-cyclic setups (pres), bugfix for correct identification of indices of extreme values in case of non-cyclic boundary conditions
Property svn:keywords set to `Id`
File size: 32.5 KB

Rev	Line
[1682]	1	!> @file pres.f90
[4649]	2	!--------------------------------------------------------------------------------------------------!
[2696]	3	! This file is part of the PALM model system.
[1036]	4	!
[4649]	5	! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
	6	! Public License as published by the Free Software Foundation, either version 3 of the License, or
	7	! (at your option) any later version.
[1036]	8	!
[4649]	9	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
	10	! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
	11	! Public License for more details.
[1036]	12	!
[4649]	13	! You should have received a copy of the GNU General Public License along with PALM. If not, see
	14	! <http://www.gnu.org/licenses/>.
[1036]	15	!
[4828]	16	! Copyright 1997-2021 Leibniz Universitaet Hannover
[4649]	17	!--------------------------------------------------------------------------------------------------!
[1036]	18	!
[4649]	19	!
[484]	20	! Current revisions:
[4649]	21	! -----------------
	22	!
	23	!
[1919]	24	! Former revisions:
	25	! -----------------
[3016]	26	! $Id: pres.f90 4855 2021-01-25 12:30:54Z raasch $
[4855]	27	! bugfix: mean w removal not applied to ghost points of the total domain in case of non-cyclic
	28	! setups
	29	!
	30	! 4828 2021-01-05 11:21:41Z Giersch
[4740]	31	! Bugfix in OpenMP directives - intel compiler do not allow reduction operations on array elements
	32	!
	33	! 4719 2020-10-01 11:28:33Z pavelkrc
[4719]	34	! Bugfix of previous commit (div_old discrepancy)
	35	!
	36	! 4717 2020-09-30 22:27:40Z pavelkrc
[4717]	37	! Fixes and optimizations of OpenMP parallelization, formatting of OpenMP
	38	! directives (J. Resler)
	39	!
	40	! 4651 2020-08-27 07:17:45Z raasch
[4651]	41	! preprocessor branch for ibm removed
	42	!
	43	! 4649 2020-08-25 12:11:17Z raasch
[4649]	44	! File re-formatted to follow the PALM coding standard
	45	!
	46	! 4457 2020-03-11 14:20:43Z raasch
[4457]	47	! use statement for exchange horiz added
[4649]	48	!
[4457]	49	! 4360 2020-01-07 11:25:50Z suehring
[4346]	50	! Introduction of wall_flags_total_0, which currently sets bits based on static
	51	! topography information used in wall_flags_static_0
[4649]	52	!
[4346]	53	! 4329 2019-12-10 15:46:36Z motisi
[4329]	54	! Renamed wall_flags_0 to wall_flags_static_0
[4649]	55	!
[4329]	56	! 4182 2019-08-22 15:20:23Z scharf
[4182]	57	! Corrected "Former revisions" section
[4649]	58	!
[4182]	59	! 4015 2019-06-05 13:25:35Z raasch
[4015]	60	! variable child_domain_nvn eliminated
[4649]	61	!
[4015]	62	! 3849 2019-04-01 16:35:16Z knoop
[3634]	63	! OpenACC port for SPEC
[1919]	64	!
[4182]	65	! Revision 1.1 1997/07/24 11:24:44 raasch
	66	! Initial revision
	67	!
	68	!
[4649]	69	!--------------------------------------------------------------------------------------------------!
[1]	70	! Description:
	71	! ------------
[4649]	72	!> Compute the divergence of the provisional velocity field. Solve the Poisson equation for the
	73	!> perturbation pressure. Compute the final velocities using this perturbation pressure. Compute the
	74	!> remaining divergence.
	75	!--------------------------------------------------------------------------------------------------!
[1682]	76	SUBROUTINE pres
[1]	77
[1320]	78
[4649]	79	USE arrays_3d, &
	80	ONLY: d, &
	81	ddzu, &
	82	ddzu_pres, &
	83	ddzw, &
	84	dzw, &
	85	p, &
	86	p_loc, &
	87	rho_air, &
	88	rho_air_zw, &
	89	tend, &
	90	u, &
	91	v, &
	92	w
[1320]	93
[4649]	94	USE control_parameters, &
	95	ONLY: bc_lr_cyc, &
	96	bc_ns_cyc, &
	97	bc_radiation_l, &
	98	bc_radiation_n, &
	99	bc_radiation_r, &
	100	bc_radiation_s, &
	101	child_domain, &
	102	conserve_volume_flow, &
	103	coupling_mode, &
	104	dt_3d, &
	105	gathered_size, &
	106	ibc_p_b, &
	107	ibc_p_t, &
	108	intermediate_timestep_count, &
	109	intermediate_timestep_count_max, &
	110	mg_switch_to_pe0_level, &
	111	nesting_offline, &
	112	psolver, &
	113	subdomain_size, &
	114	topography, &
	115	volume_flow, &
	116	volume_flow_area, &
	117	volume_flow_initial
[1320]	118
[4649]	119	USE cpulog, &
	120	ONLY: cpu_log, &
	121	log_point, &
	122	log_point_s
	123
	124	USE exchange_horiz_mod, &
[4457]	125	ONLY: exchange_horiz
	126
[4649]	127	USE grid_variables, &
	128	ONLY: ddx, &
	129	ddy
[1320]	130
[4649]	131	USE indices, &
	132	ONLY: nbgp, &
	133	ngp_2dh_outer, &
	134	nx, &
	135	nxl, &
	136	nxlg, &
	137	nxl_mg, &
	138	nxr, &
	139	nxrg, &
	140	nxr_mg, &
	141	ny, &
	142	nys, &
	143	nysg, &
	144	nys_mg, &
	145	nyn, &
	146	nyng, &
	147	nyn_mg, &
	148	nzb, &
	149	nzt, &
	150	nzt_mg, &
[4346]	151	wall_flags_total_0
[1320]	152
	153	USE kinds
	154
[1]	155	USE pegrid
	156
[4649]	157	USE pmc_interface, &
	158	ONLY: nesting_mode
	159
	160	USE poisfft_mod, &
[1320]	161	ONLY: poisfft
	162
[1575]	163	USE poismg_mod
	164
[2696]	165	USE poismg_noopt_mod
	166
[4649]	167	USE statistics, &
	168	ONLY: statistic_regions, &
	169	sums_divnew_l, &
	170	sums_divold_l, &
	171	weight_pres, &
[1320]	172	weight_substep
	173
[4649]	174	USE surface_mod, &
[2232]	175	ONLY : bc_h
	176
[1]	177	IMPLICIT NONE
	178
[4649]	179	INTEGER(iwp) :: i !<
	180	INTEGER(iwp) :: j !<
	181	INTEGER(iwp) :: k !<
	182	INTEGER(iwp) :: m !<
[1]	183
[4649]	184	REAL(wp) :: ddt_3d !<
	185	REAL(wp) :: d_weight_pres !<
[4719]	186	REAL(wp) :: localsum !<
[4649]	187	REAL(wp) :: threadsum !<
	188	REAL(wp) :: weight_pres_l !<
	189	REAL(wp) :: weight_substep_l !<
[1]	190
[1762]	191	REAL(wp), DIMENSION(1:3) :: volume_flow_l !<
	192	REAL(wp), DIMENSION(1:3) :: volume_flow_offset !<
[1682]	193	REAL(wp), DIMENSION(1:nzt) :: w_l !<
	194	REAL(wp), DIMENSION(1:nzt) :: w_l_l !<
[1]	195
	196
	197	CALL cpu_log( log_point(8), 'pres', 'start' )
	198
[1918]	199	!
	200	!-- Calculate quantities to be used locally
[1342]	201	ddt_3d = 1.0_wp / dt_3d
[1918]	202	IF ( intermediate_timestep_count == 0 ) THEN
	203	!
	204	!-- If pres is called before initial time step
[1929]	205	weight_pres_l = 1.0_wp
	206	d_weight_pres = 1.0_wp
	207	weight_substep_l = 1.0_wp
[1918]	208	ELSE
[1929]	209	weight_pres_l = weight_pres(intermediate_timestep_count)
	210	d_weight_pres = 1.0_wp / weight_pres(intermediate_timestep_count)
	211	weight_substep_l = weight_substep(intermediate_timestep_count)
[1918]	212	ENDIF
[85]	213
[1]	214	!
[707]	215	!-- Multigrid method expects array d to have one ghost layer.
[4649]	216	!--
[1575]	217	IF ( psolver(1:9) == 'multigrid' ) THEN
[4649]	218
[1]	219	DEALLOCATE( d )
[4649]	220	ALLOCATE( d(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
[707]	221
	222	!
[4649]	223	!-- Since p is later used to hold the weighted average of the substeps, it cannot be used in the
	224	!-- iterative solver. Therefore, its initial value is stored on p_loc, which is then iteratively
	225	!-- advanced in every substep.
[1762]	226	IF ( intermediate_timestep_count <= 1 ) THEN
[707]	227	DO i = nxl-1, nxr+1
	228	DO j = nys-1, nyn+1
	229	DO k = nzb, nzt+1
	230	p_loc(k,j,i) = p(k,j,i)
	231	ENDDO
	232	ENDDO
	233	ENDDO
[667]	234	ENDIF
[4649]	235
[1918]	236	ELSEIF ( psolver == 'sor' .AND. intermediate_timestep_count <= 1 ) THEN
[707]	237
	238	!
[4649]	239	!-- Since p is later used to hold the weighted average of the substeps, it cannot be used in the
	240	!-- iterative solver. Therefore, its initial value is stored on p_loc, which is then iteratively
	241	!-- advanced in every substep.
[707]	242	p_loc = p
	243
[1]	244	ENDIF
	245
	246	!
[4649]	247	!-- Conserve the volume flow at the outflow in case of non-cyclic lateral boundary conditions
	248	!-- WARNING: so far, this conservation does not work at the left/south boundary if the topography at
	249	!-- the inflow differs from that at the outflow! For this case, volume_flow_area needs
	250	!-- adjustment!
[106]	251	!
	252	!-- Left/right
[4649]	253	IF ( conserve_volume_flow .AND. ( bc_radiation_l .OR. bc_radiation_r ) ) THEN
[680]	254
[1342]	255	volume_flow(1) = 0.0_wp
	256	volume_flow_l(1) = 0.0_wp
[106]	257
[3182]	258	IF ( bc_radiation_l ) THEN
[106]	259	i = 0
[3182]	260	ELSEIF ( bc_radiation_r ) THEN
[106]	261	i = nx+1
	262	ENDIF
	263
	264	DO j = nys, nyn
	265	!
	266	!-- Sum up the volume flow through the south/north boundary
[2232]	267	DO k = nzb+1, nzt
[4649]	268	volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k) &
	269	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
[106]	270	ENDDO
	271	ENDDO
	272
[4649]	273	#if defined( __parallel )
[680]	274	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
[4649]	275	CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 1, MPI_REAL, MPI_SUM, comm1dy, ierr )
[106]	276	#else
[4649]	277	volume_flow = volume_flow_l
[106]	278	#endif
[4649]	279	volume_flow_offset(1) = ( volume_flow_initial(1) - volume_flow(1) ) / volume_flow_area(1)
[106]	280
[667]	281	DO j = nysg, nyng
[2232]	282	DO k = nzb+1, nzt
[4649]	283	u(k,j,i) = u(k,j,i) + volume_flow_offset(1) &
	284	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
[106]	285	ENDDO
	286	ENDDO
	287
	288	ENDIF
	289
	290	!
	291	!-- South/north
[3182]	292	IF ( conserve_volume_flow .AND. ( bc_radiation_n .OR. bc_radiation_s ) ) THEN
[106]	293
[1342]	294	volume_flow(2) = 0.0_wp
	295	volume_flow_l(2) = 0.0_wp
[75]	296
[3182]	297	IF ( bc_radiation_s ) THEN
[106]	298	j = 0
[3182]	299	ELSEIF ( bc_radiation_n ) THEN
[75]	300	j = ny+1
[106]	301	ENDIF
	302
	303	DO i = nxl, nxr
[75]	304	!
[106]	305	!-- Sum up the volume flow through the south/north boundary
[2232]	306	DO k = nzb+1, nzt
[4649]	307	volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k) &
	308	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
[75]	309	ENDDO
[106]	310	ENDDO
	311
[4649]	312	#if defined( __parallel )
[680]	313	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
[4649]	314	CALL MPI_ALLREDUCE( volume_flow_l(2), volume_flow(2), 1, MPI_REAL, MPI_SUM, comm1dx, ierr )
[75]	315	#else
[4649]	316	volume_flow = volume_flow_l
[75]	317	#endif
[4649]	318	volume_flow_offset(2) = ( volume_flow_initial(2) - volume_flow(2) ) / volume_flow_area(2)
[75]	319
[667]	320	DO i = nxlg, nxrg
[2232]	321	DO k = nzb+1, nzt
[4649]	322	v(k,j,i) = v(k,j,i) + volume_flow_offset(2) &
	323	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
[75]	324	ENDDO
[106]	325	ENDDO
[75]	326
	327	ENDIF
	328
[76]	329	!
[4015]	330	!-- Remove mean vertical velocity in case that Neumann conditions are used both at bottom and top
[4649]	331	!-- boundary. With Neumann conditions at both vertical boundaries, the solver cannot remove mean
	332	!-- vertical velocities. They should be removed, because incompressibility requires that the
	333	!-- vertical gradient of vertical velocity is zero. Since w=0 at the solid surface, it must be zero
	334	!-- everywhere.
[4015]	335	!-- This must not be done in case of a 3d-nesting child domain, because a mean vertical velocity
	336	!-- can physically exist in such a domain.
	337	!-- Also in case of offline nesting, mean vertical velocities may exist (and must not be removed),
	338	!-- caused by horizontal divergence/convergence of the large scale flow that is prescribed at the
	339	!-- side boundaries.
[4649]	340	!-- The removal cannot be done before the first initial time step because ngp_2dh_outer is not yet
	341	!-- known then.
[4015]	342	IF ( ibc_p_b == 1 .AND. ibc_p_t == 1 .AND. .NOT. nesting_offline &
	343	.AND. .NOT. ( child_domain .AND. nesting_mode /= 'vertical' ) &
[4649]	344	.AND. intermediate_timestep_count /= 0 ) THEN
[1918]	345	w_l = 0.0_wp; w_l_l = 0.0_wp
	346	DO i = nxl, nxr
	347	DO j = nys, nyn
[2232]	348	DO k = nzb+1, nzt
[4015]	349	w_l_l(k) = w_l_l(k) + w(k,j,i) &
[4649]	350	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
[76]	351	ENDDO
	352	ENDDO
[1918]	353	ENDDO
[4649]	354	#if defined( __parallel )
[1918]	355	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[4015]	356	CALL MPI_ALLREDUCE( w_l_l(1), w_l(1), nzt, MPI_REAL, MPI_SUM, comm2d, ierr )
[76]	357	#else
[1918]	358	w_l = w_l_l
[76]	359	#endif
[1918]	360	DO k = 1, nzt
	361	w_l(k) = w_l(k) / ngp_2dh_outer(k,0)
	362	ENDDO
[4855]	363	DO i = nxl, nxr
	364	DO j = nys, nyn
[2232]	365	DO k = nzb+1, nzt
[4015]	366	w(k,j,i) = w(k,j,i) - w_l(k) &
[4649]	367	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
[76]	368	ENDDO
	369	ENDDO
[1918]	370	ENDDO
[4855]	371	!
	372	!-- Instead of running the above loop over ghost points, they are set via exchange_horiz,
	373	!-- in order to correctly consider non-cyclic boundary conditions, where ghost boundaries
	374	!-- of the total domain must not be set. Otherwise w may continuously increase/decrease
	375	!-- at these points.
	376	CALL exchange_horiz( w, nbgp )
[76]	377	ENDIF
[75]	378
	379	!
[1]	380	!-- Compute the divergence of the provisional velocity field.
	381	CALL cpu_log( log_point_s(1), 'divergence', 'start' )
	382
[1575]	383	IF ( psolver(1:9) == 'multigrid' ) THEN
[4717]	384	!$OMP PARALLEL DO PRIVATE (i,j,k) SCHEDULE( STATIC )
[1]	385	DO i = nxl-1, nxr+1
	386	DO j = nys-1, nyn+1
	387	DO k = nzb, nzt+1
[1342]	388	d(k,j,i) = 0.0_wp
[1]	389	ENDDO
	390	ENDDO
	391	ENDDO
	392	ELSE
[4717]	393	!$OMP PARALLEL DO PRIVATE (i,j,k) SCHEDULE( STATIC )
[3634]	394	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
	395	!$ACC PRESENT(d)
[1003]	396	DO i = nxl, nxr
	397	DO j = nys, nyn
	398	DO k = nzb+1, nzt
[1342]	399	d(k,j,i) = 0.0_wp
[1]	400	ENDDO
	401	ENDDO
	402	ENDDO
	403	ENDIF
	404
[4719]	405	localsum = 0.0_wp
	406	threadsum = 0.0_wp
	407
	408	!$OMP PARALLEL PRIVATE (i,j,k)
	409	!$OMP DO SCHEDULE( STATIC )
[3634]	410	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
[4346]	411	!$ACC PRESENT(u, v, w, rho_air, rho_air_zw, ddzw, wall_flags_total_0) &
[3634]	412	!$ACC PRESENT(d)
[1111]	413	DO i = nxl, nxr
	414	DO j = nys, nyn
	415	DO k = 1, nzt
[4649]	416	d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx + &
	417	( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy + &
	418	( w(k,j,i) * rho_air_zw(k) - w(k-1,j,i) * rho_air_zw(k-1) ) &
	419	* ddzw(k) ) * ddt_3d * d_weight_pres &
	420	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
[1]	421	ENDDO
	422	ENDDO
[1111]	423	ENDDO
[4719]	424	!$OMP END PARALLEL
[1]	425
	426	!
[4649]	427	!-- Compute possible PE-sum of divergences for flow_statistics. Carry out computation only at last
	428	!-- Runge-Kutta substep.
	429	IF ( intermediate_timestep_count == intermediate_timestep_count_max .OR. &
[1918]	430	intermediate_timestep_count == 0 ) THEN
[4719]	431	!$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
	432	!$OMP DO SCHEDULE( STATIC )
[3634]	433	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
	434	!$ACC REDUCTION(+:threadsum) COPY(threadsum) &
	435	!$ACC PRESENT(d)
[1908]	436	DO i = nxl, nxr
	437	DO j = nys, nyn
	438	DO k = nzb+1, nzt
	439	threadsum = threadsum + ABS( d(k,j,i) )
	440	ENDDO
[1]	441	ENDDO
	442	ENDDO
[4719]	443	localsum = localsum + threadsum * dt_3d * weight_pres_l
	444	!$OMP END PARALLEL
[1908]	445	ENDIF
[1]	446
	447	!
[4649]	448	!-- For completeness, set the divergence sum of all statistic regions to those of the total domain
	449	IF ( intermediate_timestep_count == intermediate_timestep_count_max .OR. &
[1918]	450	intermediate_timestep_count == 0 ) THEN
[4719]	451	sums_divold_l(0:statistic_regions) = localsum
[1918]	452	ENDIF
[1]	453
	454	CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
	455
	456	!
	457	!-- Compute the pressure perturbation solving the Poisson equation
[1575]	458	IF ( psolver == 'poisfft' ) THEN
[1]	459
	460	!
	461	!-- Solve Poisson equation via FFT and solution of tridiagonal matrices
[1575]	462	CALL poisfft( d )
[1212]	463
[1]	464	!
[4649]	465	!-- Store computed perturbation pressure and set boundary condition in z-direction
[4717]	466	!$OMP PARALLEL DO PRIVATE (i,j,k) SCHEDULE( STATIC )
[3634]	467	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
	468	!$ACC PRESENT(d, tend)
[1]	469	DO i = nxl, nxr
	470	DO j = nys, nyn
	471	DO k = nzb+1, nzt
	472	tend(k,j,i) = d(k,j,i)
	473	ENDDO
	474	ENDDO
	475	ENDDO
	476
	477	!
	478	!-- Bottom boundary:
[4649]	479	!-- This condition is only required for internal output. The pressure gradient
	480	!-- (dp(nzb+1)-dp(nzb))/dz is not used anywhere else.
[1]	481	IF ( ibc_p_b == 1 ) THEN
	482	!
[4649]	483	!-- Neumann (dp/dz = 0). Using surface data type, first for non-natural surfaces, then for
	484	!-- natural and urban surfaces
[2232]	485	!-- Upward facing
[4717]	486	!$OMP PARALLEL DO PRIVATE( m, i, j, k ) SCHEDULE( STATIC )
[3634]	487	!$ACC PARALLEL LOOP PRIVATE(i, j, k) &
	488	!$ACC PRESENT(bc_h, tend)
[2232]	489	DO m = 1, bc_h(0)%ns
[2298]	490	i = bc_h(0)%i(m)
	491	j = bc_h(0)%j(m)
	492	k = bc_h(0)%k(m)
[2232]	493	tend(k-1,j,i) = tend(k,j,i)
[1]	494	ENDDO
[2232]	495	!
	496	!-- Downward facing
[4717]	497	!$OMP PARALLEL DO PRIVATE( m, i, j, k ) SCHEDULE( STATIC )
[3634]	498	!$ACC PARALLEL LOOP PRIVATE(i, j, k) &
	499	!$ACC PRESENT(bc_h, tend)
[2232]	500	DO m = 1, bc_h(1)%ns
[2298]	501	i = bc_h(1)%i(m)
	502	j = bc_h(1)%j(m)
	503	k = bc_h(1)%k(m)
[2232]	504	tend(k+1,j,i) = tend(k,j,i)
	505	ENDDO
[1]	506
	507	ELSE
	508	!
[4649]	509	!-- Dirichlet. Using surface data type, first for non-natural surfaces, then for natural and
	510	!-- urban surfaces
[2232]	511	!-- Upward facing
[4717]	512	!$OMP PARALLEL DO PRIVATE( m, i, j, k ) SCHEDULE( STATIC )
[2232]	513	DO m = 1, bc_h(0)%ns
[2298]	514	i = bc_h(0)%i(m)
	515	j = bc_h(0)%j(m)
	516	k = bc_h(0)%k(m)
[2232]	517	tend(k-1,j,i) = 0.0_wp
[1]	518	ENDDO
[2232]	519	!
	520	!-- Downward facing
[4717]	521	!$OMP PARALLEL DO PRIVATE( m, i, j, k ) SCHEDULE( STATIC )
[2232]	522	DO m = 1, bc_h(1)%ns
[2298]	523	i = bc_h(1)%i(m)
	524	j = bc_h(1)%j(m)
	525	k = bc_h(1)%k(m)
[2232]	526	tend(k+1,j,i) = 0.0_wp
	527	ENDDO
[1]	528
	529	ENDIF
	530
	531	!
	532	!-- Top boundary
	533	IF ( ibc_p_t == 1 ) THEN
	534	!
	535	!-- Neumann
[4717]	536	!$OMP PARALLEL DO PRIVATE (i,j) SCHEDULE( STATIC )
[667]	537	DO i = nxlg, nxrg
	538	DO j = nysg, nyng
[1]	539	tend(nzt+1,j,i) = tend(nzt,j,i)
	540	ENDDO
	541	ENDDO
	542
	543	ELSE
	544	!
	545	!-- Dirichlet
[4717]	546	!$OMP PARALLEL DO PRIVATE (i,j) SCHEDULE( STATIC )
[3634]	547	!$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
	548	!$ACC PRESENT(tend)
[667]	549	DO i = nxlg, nxrg
	550	DO j = nysg, nyng
[1342]	551	tend(nzt+1,j,i) = 0.0_wp
[1]	552	ENDDO
	553	ENDDO
	554
	555	ENDIF
	556
	557	!
	558	!-- Exchange boundaries for p
[667]	559	CALL exchange_horiz( tend, nbgp )
[4649]	560
[1]	561	ELSEIF ( psolver == 'sor' ) THEN
	562
	563	!
[4649]	564	!-- Solve Poisson equation for perturbation pressure using SOR-Red/Black scheme
[707]	565	CALL sor( d, ddzu_pres, ddzw, p_loc )
	566	tend = p_loc
[1]	567
[1575]	568	ELSEIF ( psolver(1:9) == 'multigrid' ) THEN
[1]	569
	570	!
[4649]	571	!-- Solve Poisson equation for perturbation pressure using Multigrid scheme, array tend is used
	572	!-- to store the residuals.
[778]	573
[4649]	574	!-- If the number of grid points of the gathered grid, which is collected on PE0, is larger than
	575	!-- the number of grid points of an PE, than array tend will be enlarged.
[778]	576	IF ( gathered_size > subdomain_size ) THEN
	577	DEALLOCATE( tend )
[4649]	578	ALLOCATE( tend(nzb:nzt_mg(mg_switch_to_pe0_level)+1,nys_mg( &
	579	mg_switch_to_pe0_level)-1:nyn_mg(mg_switch_to_pe0_level)+1, &
	580	nxl_mg(mg_switch_to_pe0_level)-1:nxr_mg(mg_switch_to_pe0_level)+1) )
[778]	581	ENDIF
	582
[1575]	583	IF ( psolver == 'multigrid' ) THEN
	584	CALL poismg( tend )
	585	ELSE
[1931]	586	CALL poismg_noopt( tend )
[1575]	587	ENDIF
[707]	588
[778]	589	IF ( gathered_size > subdomain_size ) THEN
	590	DEALLOCATE( tend )
	591	ALLOCATE( tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	592	ENDIF
	593
[1]	594	!
[4649]	595	!-- Restore perturbation pressure on tend because this array is used further below to correct the
	596	!-- velocity fields
[707]	597	DO i = nxl-1, nxr+1
	598	DO j = nys-1, nyn+1
	599	DO k = nzb, nzt+1
	600	tend(k,j,i) = p_loc(k,j,i)
	601	ENDDO
	602	ENDDO
	603	ENDDO
[667]	604
[1]	605	ENDIF
	606
	607	!
[707]	608	!-- Store perturbation pressure on array p, used for pressure data output.
	609	!-- Ghost layers are added in the output routines (except sor-method: see below)
[1918]	610	IF ( intermediate_timestep_count <= 1 ) THEN
[4717]	611	!$OMP PARALLEL DO PRIVATE (i,j,k) SCHEDULE( STATIC )
[3634]	612	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
	613	!$ACC PRESENT(p, tend)
[707]	614	DO i = nxl-1, nxr+1
	615	DO j = nys-1, nyn+1
	616	DO k = nzb, nzt+1
[4649]	617	p(k,j,i) = tend(k,j,i) * weight_substep_l
[673]	618	ENDDO
[1]	619	ENDDO
[707]	620	ENDDO
	621
[1762]	622	ELSEIF ( intermediate_timestep_count > 1 ) THEN
[4717]	623	!$OMP PARALLEL DO PRIVATE (i,j,k) SCHEDULE( STATIC )
[3634]	624	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
	625	!$ACC PRESENT(p, tend)
[707]	626	DO i = nxl-1, nxr+1
	627	DO j = nys-1, nyn+1
	628	DO k = nzb, nzt+1
[4649]	629	p(k,j,i) = p(k,j,i) + tend(k,j,i) * weight_substep_l
[673]	630	ENDDO
	631	ENDDO
[707]	632	ENDDO
	633
	634	ENDIF
[4649]	635
[707]	636	!
	637	!-- SOR-method needs ghost layers for the next timestep
	638	IF ( psolver == 'sor' ) CALL exchange_horiz( p, nbgp )
[682]	639
[1]	640	!
[4649]	641	!-- Correction of the provisional velocities with the current perturbation pressure just computed
[709]	642	IF ( conserve_volume_flow .AND. ( bc_lr_cyc .OR. bc_ns_cyc ) ) THEN
[1342]	643	volume_flow_l(1) = 0.0_wp
	644	volume_flow_l(2) = 0.0_wp
[1]	645	ENDIF
[3347]	646	!
[4649]	647	!-- Add pressure gradients to the velocity components. Note, the loops are running over the entire
	648	!-- model domain, even though, in case of non-cyclic boundaries u- and v-component are not
	649	!-- prognostic at i=0 and j=0, respectiveley. However, in case of Dirichlet boundary conditions for
	650	!-- the velocities, zero-gradient conditions for the pressure are set, so that no modification is
	651	!-- imposed at the boundaries.
[4717]	652	!$OMP PARALLEL DO PRIVATE (i,j,k) SCHEDULE( STATIC )
[3634]	653	!$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k) &
[4346]	654	!$ACC PRESENT(u, v, w, tend, ddzu, wall_flags_total_0)
[4649]	655	DO i = nxl, nxr
[1]	656	DO j = nys, nyn
[2118]	657
[2232]	658	DO k = nzb+1, nzt
[4649]	659	w(k,j,i) = w(k,j,i) - dt_3d * ( tend(k+1,j,i) - tend(k,j,i) ) * ddzu(k+1) &
	660	* weight_pres_l &
	661	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
[1]	662	ENDDO
[2118]	663
[2232]	664	DO k = nzb+1, nzt
[4649]	665	u(k,j,i) = u(k,j,i) - dt_3d * ( tend(k,j,i) - tend(k,j,i-1) ) * ddx * weight_pres_l &
	666	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
[1]	667	ENDDO
[2118]	668
[2232]	669	DO k = nzb+1, nzt
[4649]	670	v(k,j,i) = v(k,j,i) - dt_3d * ( tend(k,j,i) - tend(k,j-1,i) ) * ddy * weight_pres_l &
	671	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
	672	ENDDO
[1]	673
	674	ENDDO
	675	ENDDO
[1113]	676
	677	!
[4649]	678	!-- The vertical velocity is not set to zero at nzt + 1 for nested domains. Instead it is set to the
	679	!-- values of nzt (see routine vnest_boundary_conds or pmci_interp_tril_t) BEFORE calling the
	680	!-- pressure solver. To avoid jumps while plotting profiles, w at the top has to be set to the
	681	!-- values in height nzt after above modifications. Hint: w level nzt+1 does not impact results.
	682	IF ( child_domain .OR. coupling_mode == 'vnested_fine' ) THEN
[3012]	683	w(nzt+1,:,:) = w(nzt,:,:)
	684	ENDIF
	685
	686	!
[1113]	687	!-- Sum up the volume flow through the right and north boundary
[4649]	688	IF ( conserve_volume_flow .AND. bc_lr_cyc .AND. bc_ns_cyc .AND. nxr == nx ) THEN
[4740]	689	threadsum = 0.0_wp
	690	!
	691	!-- Summation of the volume flow is done on threadsum rather than on volumen_flow itself.
	692	!-- This is because intel compiler, when compiled with openmp, do not allow reduction
	693	!-- operation on array elements.
	694	!$OMP PARALLEL DO PRIVATE (j,k) REDUCTION (+:threadsum)
[1113]	695	DO j = nys, nyn
[2232]	696	DO k = nzb+1, nzt
[4740]	697	threadsum = threadsum + u(k,j,nxr) * dzw(k) &
[4649]	698	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,nxr), 1 ) )
[1113]	699	ENDDO
	700	ENDDO
[4740]	701	volume_flow_l(1) = threadsum
[1113]	702
	703	ENDIF
	704
[4649]	705	IF ( conserve_volume_flow .AND. bc_ns_cyc .AND. bc_lr_cyc .AND. nyn == ny ) THEN
[4740]	706	threadsum = 0.0_wp
	707	!$OMP PARALLEL DO PRIVATE (j,k) REDUCTION (+:threadsum)
[1113]	708	DO i = nxl, nxr
[2232]	709	DO k = nzb+1, nzt
[4740]	710	threadsum = threadsum + v(k,nyn,i) * dzw(k) &
[4649]	711	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,nyn,i), 2 ) )
[1113]	712	ENDDO
	713	ENDDO
[4740]	714	volume_flow_l(2) = threadsum
[1113]	715
	716	ENDIF
[4649]	717
[1]	718	!
	719	!-- Conserve the volume flow
[707]	720	IF ( conserve_volume_flow .AND. ( bc_lr_cyc .AND. bc_ns_cyc ) ) THEN
[1]	721
[4649]	722	#if defined( __parallel )
[622]	723	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[4649]	724	CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 2, MPI_REAL, MPI_SUM, comm2d, ierr )
[1]	725	#else
[4649]	726	volume_flow = volume_flow_l
	727	#endif
[1]	728
[4649]	729	volume_flow_offset(1:2) = ( volume_flow_initial(1:2) - volume_flow(1:2) ) &
	730	/ volume_flow_area(1:2)
[1]	731
[4717]	732	!$OMP PARALLEL DO PRIVATE (i,j,k) SCHEDULE( STATIC )
[1]	733	DO i = nxl, nxr
	734	DO j = nys, nyn
[2232]	735	DO k = nzb+1, nzt
[4649]	736	u(k,j,i) = u(k,j,i) + volume_flow_offset(1) &
	737	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
[719]	738	ENDDO
[2232]	739	DO k = nzb+1, nzt
[4649]	740	v(k,j,i) = v(k,j,i) + volume_flow_offset(2) &
	741	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
[667]	742	ENDDO
[1]	743	ENDDO
	744	ENDDO
[667]	745
[1]	746	ENDIF
	747
	748	!
	749	!-- Exchange of boundaries for the velocities
[667]	750	CALL exchange_horiz( u, nbgp )
	751	CALL exchange_horiz( v, nbgp )
	752	CALL exchange_horiz( w, nbgp )
[1]	753
	754	!
[4649]	755	!-- Compute the divergence of the corrected velocity field.
	756	!-- A possible PE-sum is computed in flow_statistics. Carry out computation only at last
	757	!-- Runge-Kutta step.
	758	IF ( intermediate_timestep_count == intermediate_timestep_count_max .OR. &
[1918]	759	intermediate_timestep_count == 0 ) THEN
[1908]	760	CALL cpu_log( log_point_s(1), 'divergence', 'start' )
	761	sums_divnew_l = 0.0_wp
[1]	762
	763	!
[4649]	764	!-- d must be reset to zero because it can contain nonzero values below the topography
[1908]	765	IF ( topography /= 'flat' ) d = 0.0_wp
[1]	766
[4719]	767	localsum = 0.0_wp
	768	threadsum = 0.0_wp
	769
	770	!$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
	771	!$OMP DO SCHEDULE( STATIC )
[3634]	772	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
[4346]	773	!$ACC PRESENT(u, v, w, rho_air, rho_air_zw, ddzw, wall_flags_total_0) &
[3634]	774	!$ACC PRESENT(d)
[1908]	775	DO i = nxl, nxr
	776	DO j = nys, nyn
[2232]	777	DO k = nzb+1, nzt
[2037]	778	d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx + &
	779	( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy + &
[4649]	780	( w(k,j,i) * rho_air_zw(k) - w(k-1,j,i) * rho_air_zw(k-1) ) &
	781	* ddzw(k) ) &
	782	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
[1908]	783	ENDDO
[1113]	784	ENDDO
	785	ENDDO
	786	!
[1908]	787	!-- Compute possible PE-sum of divergences for flow_statistics
[4719]	788	!$OMP DO SCHEDULE( STATIC )
[3634]	789	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
	790	!$ACC REDUCTION(+:threadsum) COPY(threadsum) &
	791	!$ACC PRESENT(d)
[1908]	792	DO i = nxl, nxr
	793	DO j = nys, nyn
	794	DO k = nzb+1, nzt
	795	threadsum = threadsum + ABS( d(k,j,i) )
	796	ENDDO
[1]	797	ENDDO
	798	ENDDO
[667]	799
[4719]	800	localsum = localsum + threadsum
	801	!$OMP END PARALLEL
	802
[1]	803	!
[4649]	804	!-- For completeness, set the divergence sum of all statistic regions to those of the total
	805	!-- domain
[4719]	806	sums_divnew_l(0:statistic_regions) = localsum
[1]	807
[1908]	808	CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
[1]	809
[1908]	810	ENDIF
	811
[1]	812	CALL cpu_log( log_point(8), 'pres', 'stop' )
	813
[4651]	814	END SUBROUTINE pres

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

Context Navigation

source: palm/trunk/SOURCE/pres.f90 @ 4899

Download in other formats: