Home

Context Navigation

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

Last change on this file since 4717 was 4717, checked in by pavelkrc, 4 years ago
Fixes and optimizations of OpenMP parallelization, formatting of OpenMP directives
Property svn:keywords set to `Id`
File size: 31.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |