Home

Context Navigation

pres.f90 @ 4879

Last change on this file since 4879 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

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-2021 Leibniz Universitaet Hannover
17	!--------------------------------------------------------------------------------------------------!
18	!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: pres.f90 4855 2021-01-25 12:30:54Z gronemeier $
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
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
34	! Bugfix of previous commit (div_old discrepancy)
35	!
36	! 4717 2020-09-30 22:27:40Z pavelkrc
37	! Fixes and optimizations of OpenMP parallelization, formatting of OpenMP
38	! directives (J. Resler)
39	!
40	! 4651 2020-08-27 07:17:45Z raasch
41	! preprocessor branch for ibm removed
42	!
43	! 4649 2020-08-25 12:11:17Z raasch
44	! File re-formatted to follow the PALM coding standard
45	!
46	! 4457 2020-03-11 14:20:43Z raasch
47	! use statement for exchange horiz added
48	!
49	! 4360 2020-01-07 11:25:50Z suehring
50	! Introduction of wall_flags_total_0, which currently sets bits based on static
51	! topography information used in wall_flags_static_0
52	!
53	! 4329 2019-12-10 15:46:36Z motisi
54	! Renamed wall_flags_0 to wall_flags_static_0
55	!
56	! 4182 2019-08-22 15:20:23Z scharf
57	! Corrected "Former revisions" section
58	!
59	! 4015 2019-06-05 13:25:35Z raasch
60	! variable child_domain_nvn eliminated
61	!
62	! 3849 2019-04-01 16:35:16Z knoop
63	! OpenACC port for SPEC
64	!
65	! Revision 1.1 1997/07/24 11:24:44 raasch
66	! Initial revision
67	!
68	!
69	!--------------------------------------------------------------------------------------------------!
70	! Description:
71	! ------------
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	!--------------------------------------------------------------------------------------------------!
76	SUBROUTINE pres
77
78
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
93
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
118
119	USE cpulog, &
120	ONLY: cpu_log, &
121	log_point, &
122	log_point_s
123
124	USE exchange_horiz_mod, &
125	ONLY: exchange_horiz
126
127	USE grid_variables, &
128	ONLY: ddx, &
129	ddy
130
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, &
151	wall_flags_total_0
152
153	USE kinds
154
155	USE pegrid
156
157	USE pmc_interface, &
158	ONLY: nesting_mode
159
160	USE poisfft_mod, &
161	ONLY: poisfft
162
163	USE poismg_mod
164
165	USE poismg_noopt_mod
166
167	USE statistics, &
168	ONLY: statistic_regions, &
169	sums_divnew_l, &
170	sums_divold_l, &
171	weight_pres, &
172	weight_substep
173
174	USE surface_mod, &
175	ONLY : bc_h
176
177	IMPLICIT NONE
178
179	INTEGER(iwp) :: i !<
180	INTEGER(iwp) :: j !<
181	INTEGER(iwp) :: k !<
182	INTEGER(iwp) :: m !<
183
184	REAL(wp) :: ddt_3d !<
185	REAL(wp) :: d_weight_pres !<
186	REAL(wp) :: localsum !<
187	REAL(wp) :: threadsum !<
188	REAL(wp) :: weight_pres_l !<
189	REAL(wp) :: weight_substep_l !<
190
191	REAL(wp), DIMENSION(1:3) :: volume_flow_l !<
192	REAL(wp), DIMENSION(1:3) :: volume_flow_offset !<
193	REAL(wp), DIMENSION(1:nzt) :: w_l !<
194	REAL(wp), DIMENSION(1:nzt) :: w_l_l !<
195
196
197	CALL cpu_log( log_point(8), 'pres', 'start' )
198
199	!
200	!-- Calculate quantities to be used locally
201	ddt_3d = 1.0_wp / dt_3d
202	IF ( intermediate_timestep_count == 0 ) THEN
203	!
204	!-- If pres is called before initial time step
205	weight_pres_l = 1.0_wp
206	d_weight_pres = 1.0_wp
207	weight_substep_l = 1.0_wp
208	ELSE
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)
212	ENDIF
213
214	!
215	!-- Multigrid method expects array d to have one ghost layer.
216	!--
217	IF ( psolver(1:9) == 'multigrid' ) THEN
218
219	DEALLOCATE( d )
220	ALLOCATE( d(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
221
222	!
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.
226	IF ( intermediate_timestep_count <= 1 ) THEN
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
234	ENDIF
235
236	ELSEIF ( psolver == 'sor' .AND. intermediate_timestep_count <= 1 ) THEN
237
238	!
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.
242	p_loc = p
243
244	ENDIF
245
246	!
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!
251	!
252	!-- Left/right
253	IF ( conserve_volume_flow .AND. ( bc_radiation_l .OR. bc_radiation_r ) ) THEN
254
255	volume_flow(1) = 0.0_wp
256	volume_flow_l(1) = 0.0_wp
257
258	IF ( bc_radiation_l ) THEN
259	i = 0
260	ELSEIF ( bc_radiation_r ) THEN
261	i = nx+1
262	ENDIF
263
264	DO j = nys, nyn
265	!
266	!-- Sum up the volume flow through the south/north boundary
267	DO k = nzb+1, nzt
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 ) )
270	ENDDO
271	ENDDO
272
273	#if defined( __parallel )
274	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
275	CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 1, MPI_REAL, MPI_SUM, comm1dy, ierr )
276	#else
277	volume_flow = volume_flow_l
278	#endif
279	volume_flow_offset(1) = ( volume_flow_initial(1) - volume_flow(1) ) / volume_flow_area(1)
280
281	DO j = nysg, nyng
282	DO k = nzb+1, nzt
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 ) )
285	ENDDO
286	ENDDO
287
288	ENDIF
289
290	!
291	!-- South/north
292	IF ( conserve_volume_flow .AND. ( bc_radiation_n .OR. bc_radiation_s ) ) THEN
293
294	volume_flow(2) = 0.0_wp
295	volume_flow_l(2) = 0.0_wp
296
297	IF ( bc_radiation_s ) THEN
298	j = 0
299	ELSEIF ( bc_radiation_n ) THEN
300	j = ny+1
301	ENDIF
302
303	DO i = nxl, nxr
304	!
305	!-- Sum up the volume flow through the south/north boundary
306	DO k = nzb+1, nzt
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 ) )
309	ENDDO
310	ENDDO
311
312	#if defined( __parallel )
313	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
314	CALL MPI_ALLREDUCE( volume_flow_l(2), volume_flow(2), 1, MPI_REAL, MPI_SUM, comm1dx, ierr )
315	#else
316	volume_flow = volume_flow_l
317	#endif
318	volume_flow_offset(2) = ( volume_flow_initial(2) - volume_flow(2) ) / volume_flow_area(2)
319
320	DO i = nxlg, nxrg
321	DO k = nzb+1, nzt
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 ) )
324	ENDDO
325	ENDDO
326
327	ENDIF
328
329	!
330	!-- Remove mean vertical velocity in case that Neumann conditions are used both at bottom and top
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.
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.
340	!-- The removal cannot be done before the first initial time step because ngp_2dh_outer is not yet
341	!-- known then.
342	IF ( ibc_p_b == 1 .AND. ibc_p_t == 1 .AND. .NOT. nesting_offline &
343	.AND. .NOT. ( child_domain .AND. nesting_mode /= 'vertical' ) &
344	.AND. intermediate_timestep_count /= 0 ) THEN
345	w_l = 0.0_wp; w_l_l = 0.0_wp
346	DO i = nxl, nxr
347	DO j = nys, nyn
348	DO k = nzb+1, nzt
349	w_l_l(k) = w_l_l(k) + w(k,j,i) &
350	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
351	ENDDO
352	ENDDO
353	ENDDO
354	#if defined( __parallel )
355	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
356	CALL MPI_ALLREDUCE( w_l_l(1), w_l(1), nzt, MPI_REAL, MPI_SUM, comm2d, ierr )
357	#else
358	w_l = w_l_l
359	#endif
360	DO k = 1, nzt
361	w_l(k) = w_l(k) / ngp_2dh_outer(k,0)
362	ENDDO
363	DO i = nxl, nxr
364	DO j = nys, nyn
365	DO k = nzb+1, nzt
366	w(k,j,i) = w(k,j,i) - w_l(k) &
367	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
368	ENDDO
369	ENDDO
370	ENDDO
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 )
377	ENDIF
378
379	!
380	!-- Compute the divergence of the provisional velocity field.
381	CALL cpu_log( log_point_s(1), 'divergence', 'start' )
382
383	IF ( psolver(1:9) == 'multigrid' ) THEN
384	!$OMP PARALLEL DO PRIVATE (i,j,k) SCHEDULE( STATIC )
385	DO i = nxl-1, nxr+1
386	DO j = nys-1, nyn+1
387	DO k = nzb, nzt+1
388	d(k,j,i) = 0.0_wp
389	ENDDO
390	ENDDO
391	ENDDO
392	ELSE
393	!$OMP PARALLEL DO PRIVATE (i,j,k) SCHEDULE( STATIC )
394	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
395	!$ACC PRESENT(d)
396	DO i = nxl, nxr
397	DO j = nys, nyn
398	DO k = nzb+1, nzt
399	d(k,j,i) = 0.0_wp
400	ENDDO
401	ENDDO
402	ENDDO
403	ENDIF
404
405	localsum = 0.0_wp
406	threadsum = 0.0_wp
407
408	!$OMP PARALLEL PRIVATE (i,j,k)
409	!$OMP DO SCHEDULE( STATIC )
410	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
411	!$ACC PRESENT(u, v, w, rho_air, rho_air_zw, ddzw, wall_flags_total_0) &
412	!$ACC PRESENT(d)
413	DO i = nxl, nxr
414	DO j = nys, nyn
415	DO k = 1, nzt
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 ) )
421	ENDDO
422	ENDDO
423	ENDDO
424	!$OMP END PARALLEL
425
426	!
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. &
430	intermediate_timestep_count == 0 ) THEN
431	!$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
432	!$OMP DO SCHEDULE( STATIC )
433	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
434	!$ACC REDUCTION(+:threadsum) COPY(threadsum) &
435	!$ACC PRESENT(d)
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
441	ENDDO
442	ENDDO
443	localsum = localsum + threadsum * dt_3d * weight_pres_l
444	!$OMP END PARALLEL
445	ENDIF
446
447	!
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. &
450	intermediate_timestep_count == 0 ) THEN
451	sums_divold_l(0:statistic_regions) = localsum
452	ENDIF
453
454	CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
455
456	!
457	!-- Compute the pressure perturbation solving the Poisson equation
458	IF ( psolver == 'poisfft' ) THEN
459
460	!
461	!-- Solve Poisson equation via FFT and solution of tridiagonal matrices
462	CALL poisfft( d )
463
464	!
465	!-- Store computed perturbation pressure and set boundary condition in z-direction
466	!$OMP PARALLEL DO PRIVATE (i,j,k) SCHEDULE( STATIC )
467	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
468	!$ACC PRESENT(d, tend)
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:
479	!-- This condition is only required for internal output. The pressure gradient
480	!-- (dp(nzb+1)-dp(nzb))/dz is not used anywhere else.
481	IF ( ibc_p_b == 1 ) THEN
482	!
483	!-- Neumann (dp/dz = 0). Using surface data type, first for non-natural surfaces, then for
484	!-- natural and urban surfaces
485	!-- Upward facing
486	!$OMP PARALLEL DO PRIVATE( m, i, j, k ) SCHEDULE( STATIC )
487	!$ACC PARALLEL LOOP PRIVATE(i, j, k) &
488	!$ACC PRESENT(bc_h, tend)
489	DO m = 1, bc_h(0)%ns
490	i = bc_h(0)%i(m)
491	j = bc_h(0)%j(m)
492	k = bc_h(0)%k(m)
493	tend(k-1,j,i) = tend(k,j,i)
494	ENDDO
495	!
496	!-- Downward facing
497	!$OMP PARALLEL DO PRIVATE( m, i, j, k ) SCHEDULE( STATIC )
498	!$ACC PARALLEL LOOP PRIVATE(i, j, k) &
499	!$ACC PRESENT(bc_h, tend)
500	DO m = 1, bc_h(1)%ns
501	i = bc_h(1)%i(m)
502	j = bc_h(1)%j(m)
503	k = bc_h(1)%k(m)
504	tend(k+1,j,i) = tend(k,j,i)
505	ENDDO
506
507	ELSE
508	!
509	!-- Dirichlet. Using surface data type, first for non-natural surfaces, then for natural and
510	!-- urban surfaces
511	!-- Upward facing
512	!$OMP PARALLEL DO PRIVATE( m, i, j, k ) SCHEDULE( STATIC )
513	DO m = 1, bc_h(0)%ns
514	i = bc_h(0)%i(m)
515	j = bc_h(0)%j(m)
516	k = bc_h(0)%k(m)
517	tend(k-1,j,i) = 0.0_wp
518	ENDDO
519	!
520	!-- Downward facing
521	!$OMP PARALLEL DO PRIVATE( m, i, j, k ) SCHEDULE( STATIC )
522	DO m = 1, bc_h(1)%ns
523	i = bc_h(1)%i(m)
524	j = bc_h(1)%j(m)
525	k = bc_h(1)%k(m)
526	tend(k+1,j,i) = 0.0_wp
527	ENDDO
528
529	ENDIF
530
531	!
532	!-- Top boundary
533	IF ( ibc_p_t == 1 ) THEN
534	!
535	!-- Neumann
536	!$OMP PARALLEL DO PRIVATE (i,j) SCHEDULE( STATIC )
537	DO i = nxlg, nxrg
538	DO j = nysg, nyng
539	tend(nzt+1,j,i) = tend(nzt,j,i)
540	ENDDO
541	ENDDO
542
543	ELSE
544	!
545	!-- Dirichlet
546	!$OMP PARALLEL DO PRIVATE (i,j) SCHEDULE( STATIC )
547	!$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
548	!$ACC PRESENT(tend)
549	DO i = nxlg, nxrg
550	DO j = nysg, nyng
551	tend(nzt+1,j,i) = 0.0_wp
552	ENDDO
553	ENDDO
554
555	ENDIF
556
557	!
558	!-- Exchange boundaries for p
559	CALL exchange_horiz( tend, nbgp )
560
561	ELSEIF ( psolver == 'sor' ) THEN
562
563	!
564	!-- Solve Poisson equation for perturbation pressure using SOR-Red/Black scheme
565	CALL sor( d, ddzu_pres, ddzw, p_loc )
566	tend = p_loc
567
568	ELSEIF ( psolver(1:9) == 'multigrid' ) THEN
569
570	!
571	!-- Solve Poisson equation for perturbation pressure using Multigrid scheme, array tend is used
572	!-- to store the residuals.
573
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.
576	IF ( gathered_size > subdomain_size ) THEN
577	DEALLOCATE( tend )
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) )
581	ENDIF
582
583	IF ( psolver == 'multigrid' ) THEN
584	CALL poismg( tend )
585	ELSE
586	CALL poismg_noopt( tend )
587	ENDIF
588
589	IF ( gathered_size > subdomain_size ) THEN
590	DEALLOCATE( tend )
591	ALLOCATE( tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
592	ENDIF
593
594	!
595	!-- Restore perturbation pressure on tend because this array is used further below to correct the
596	!-- velocity fields
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
604
605	ENDIF
606
607	!
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)
610	IF ( intermediate_timestep_count <= 1 ) THEN
611	!$OMP PARALLEL DO PRIVATE (i,j,k) SCHEDULE( STATIC )
612	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
613	!$ACC PRESENT(p, tend)
614	DO i = nxl-1, nxr+1
615	DO j = nys-1, nyn+1
616	DO k = nzb, nzt+1
617	p(k,j,i) = tend(k,j,i) * weight_substep_l
618	ENDDO
619	ENDDO
620	ENDDO
621
622	ELSEIF ( intermediate_timestep_count > 1 ) THEN
623	!$OMP PARALLEL DO PRIVATE (i,j,k) SCHEDULE( STATIC )
624	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
625	!$ACC PRESENT(p, tend)
626	DO i = nxl-1, nxr+1
627	DO j = nys-1, nyn+1
628	DO k = nzb, nzt+1
629	p(k,j,i) = p(k,j,i) + tend(k,j,i) * weight_substep_l
630	ENDDO
631	ENDDO
632	ENDDO
633
634	ENDIF
635
636	!
637	!-- SOR-method needs ghost layers for the next timestep
638	IF ( psolver == 'sor' ) CALL exchange_horiz( p, nbgp )
639
640	!
641	!-- Correction of the provisional velocities with the current perturbation pressure just computed
642	IF ( conserve_volume_flow .AND. ( bc_lr_cyc .OR. bc_ns_cyc ) ) THEN
643	volume_flow_l(1) = 0.0_wp
644	volume_flow_l(2) = 0.0_wp
645	ENDIF
646	!
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.
652	!$OMP PARALLEL DO PRIVATE (i,j,k) SCHEDULE( STATIC )
653	!$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k) &
654	!$ACC PRESENT(u, v, w, tend, ddzu, wall_flags_total_0)
655	DO i = nxl, nxr
656	DO j = nys, nyn
657
658	DO k = nzb+1, nzt
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 ) )
662	ENDDO
663
664	DO k = nzb+1, nzt
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 ) )
667	ENDDO
668
669	DO k = nzb+1, nzt
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
673
674	ENDDO
675	ENDDO
676
677	!
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
683	w(nzt+1,:,:) = w(nzt,:,:)
684	ENDIF
685
686	!
687	!-- Sum up the volume flow through the right and north boundary
688	IF ( conserve_volume_flow .AND. bc_lr_cyc .AND. bc_ns_cyc .AND. nxr == nx ) THEN
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)
695	DO j = nys, nyn
696	DO k = nzb+1, nzt
697	threadsum = threadsum + u(k,j,nxr) * dzw(k) &
698	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,nxr), 1 ) )
699	ENDDO
700	ENDDO
701	volume_flow_l(1) = threadsum
702
703	ENDIF
704
705	IF ( conserve_volume_flow .AND. bc_ns_cyc .AND. bc_lr_cyc .AND. nyn == ny ) THEN
706	threadsum = 0.0_wp
707	!$OMP PARALLEL DO PRIVATE (j,k) REDUCTION (+:threadsum)
708	DO i = nxl, nxr
709	DO k = nzb+1, nzt
710	threadsum = threadsum + v(k,nyn,i) * dzw(k) &
711	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,nyn,i), 2 ) )
712	ENDDO
713	ENDDO
714	volume_flow_l(2) = threadsum
715
716	ENDIF
717
718	!
719	!-- Conserve the volume flow
720	IF ( conserve_volume_flow .AND. ( bc_lr_cyc .AND. bc_ns_cyc ) ) THEN
721
722	#if defined( __parallel )
723	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
724	CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 2, MPI_REAL, MPI_SUM, comm2d, ierr )
725	#else
726	volume_flow = volume_flow_l
727	#endif
728
729	volume_flow_offset(1:2) = ( volume_flow_initial(1:2) - volume_flow(1:2) ) &
730	/ volume_flow_area(1:2)
731
732	!$OMP PARALLEL DO PRIVATE (i,j,k) SCHEDULE( STATIC )
733	DO i = nxl, nxr
734	DO j = nys, nyn
735	DO k = nzb+1, nzt
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 ) )
738	ENDDO
739	DO k = nzb+1, nzt
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 ) )
742	ENDDO
743	ENDDO
744	ENDDO
745
746	ENDIF
747
748	!
749	!-- Exchange of boundaries for the velocities
750	CALL exchange_horiz( u, nbgp )
751	CALL exchange_horiz( v, nbgp )
752	CALL exchange_horiz( w, nbgp )
753
754	!
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. &
759	intermediate_timestep_count == 0 ) THEN
760	CALL cpu_log( log_point_s(1), 'divergence', 'start' )
761	sums_divnew_l = 0.0_wp
762
763	!
764	!-- d must be reset to zero because it can contain nonzero values below the topography
765	IF ( topography /= 'flat' ) d = 0.0_wp
766
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 )
772	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
773	!$ACC PRESENT(u, v, w, rho_air, rho_air_zw, ddzw, wall_flags_total_0) &
774	!$ACC PRESENT(d)
775	DO i = nxl, nxr
776	DO j = nys, nyn
777	DO k = nzb+1, nzt
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 + &
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 ) )
783	ENDDO
784	ENDDO
785	ENDDO
786	!
787	!-- Compute possible PE-sum of divergences for flow_statistics
788	!$OMP DO SCHEDULE( STATIC )
789	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
790	!$ACC REDUCTION(+:threadsum) COPY(threadsum) &
791	!$ACC PRESENT(d)
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
797	ENDDO
798	ENDDO
799
800	localsum = localsum + threadsum
801	!$OMP END PARALLEL
802
803	!
804	!-- For completeness, set the divergence sum of all statistic regions to those of the total
805	!-- domain
806	sums_divnew_l(0:statistic_regions) = localsum
807
808	CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
809
810	ENDIF
811
812	CALL cpu_log( log_point(8), 'pres', 'stop' )
813
814	END SUBROUTINE pres

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

Context Navigation

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

Download in other formats: