Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |