Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |