Home

Context Navigation

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

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

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 4719 2020-10-01 11:28:33Z raasch $
27	! Bugfix of previous commit (div_old discrepancy)
28	!
29	! 4717 2020-09-30 22:27:40Z pavelkrc
30	! Fixes and optimizations of OpenMP parallelization, formatting of OpenMP
31	! directives (J. Resler)
32	!
33	! 4651 2020-08-27 07:17:45Z raasch
34	! preprocessor branch for ibm removed
35	!
36	! 4649 2020-08-25 12:11:17Z raasch
37	! File re-formatted to follow the PALM coding standard
38	!
39	! 4457 2020-03-11 14:20:43Z raasch
40	! use statement for exchange horiz added
41	!
42	! 4360 2020-01-07 11:25:50Z suehring
43	! Introduction of wall_flags_total_0, which currently sets bits based on static
44	! topography information used in wall_flags_static_0
45	!
46	! 4329 2019-12-10 15:46:36Z motisi
47	! Renamed wall_flags_0 to wall_flags_static_0
48	!
49	! 4182 2019-08-22 15:20:23Z scharf
50	! Corrected "Former revisions" section
51	!
52	! 4015 2019-06-05 13:25:35Z raasch
53	! variable child_domain_nvn eliminated
54	!
55	! 3849 2019-04-01 16:35:16Z knoop
56	! OpenACC port for SPEC
57	!
58	! Revision 1.1 1997/07/24 11:24:44 raasch
59	! Initial revision
60	!
61	!
62	!--------------------------------------------------------------------------------------------------!
63	! Description:
64	! ------------
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	!--------------------------------------------------------------------------------------------------!
69	SUBROUTINE pres
70
71
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
86
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
111
112	USE cpulog, &
113	ONLY: cpu_log, &
114	log_point, &
115	log_point_s
116
117	USE exchange_horiz_mod, &
118	ONLY: exchange_horiz
119
120	USE grid_variables, &
121	ONLY: ddx, &
122	ddy
123
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, &
144	wall_flags_total_0
145
146	USE kinds
147
148	USE pegrid
149
150	USE pmc_interface, &
151	ONLY: nesting_mode
152
153	USE poisfft_mod, &
154	ONLY: poisfft
155
156	USE poismg_mod
157
158	USE poismg_noopt_mod
159
160	USE statistics, &
161	ONLY: statistic_regions, &
162	sums_divnew_l, &
163	sums_divold_l, &
164	weight_pres, &
165	weight_substep
166
167	USE surface_mod, &
168	ONLY : bc_h
169
170	IMPLICIT NONE
171
172	INTEGER(iwp) :: i !<
173	INTEGER(iwp) :: j !<
174	INTEGER(iwp) :: k !<
175	INTEGER(iwp) :: m !<
176
177	REAL(wp) :: ddt_3d !<
178	REAL(wp) :: d_weight_pres !<
179	REAL(wp) :: localsum !<
180	REAL(wp) :: threadsum !<
181	REAL(wp) :: weight_pres_l !<
182	REAL(wp) :: weight_substep_l !<
183
184	REAL(wp), DIMENSION(1:3) :: volume_flow_l !<
185	REAL(wp), DIMENSION(1:3) :: volume_flow_offset !<
186	REAL(wp), DIMENSION(1:nzt) :: w_l !<
187	REAL(wp), DIMENSION(1:nzt) :: w_l_l !<
188
189
190	CALL cpu_log( log_point(8), 'pres', 'start' )
191
192	!
193	!-- Calculate quantities to be used locally
194	ddt_3d = 1.0_wp / dt_3d
195	IF ( intermediate_timestep_count == 0 ) THEN
196	!
197	!-- If pres is called before initial time step
198	weight_pres_l = 1.0_wp
199	d_weight_pres = 1.0_wp
200	weight_substep_l = 1.0_wp
201	ELSE
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)
205	ENDIF
206
207	!
208	!-- Multigrid method expects array d to have one ghost layer.
209	!--
210	IF ( psolver(1:9) == 'multigrid' ) THEN
211
212	DEALLOCATE( d )
213	ALLOCATE( d(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
214
215	!
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.
219	IF ( intermediate_timestep_count <= 1 ) THEN
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
227	ENDIF
228
229	ELSEIF ( psolver == 'sor' .AND. intermediate_timestep_count <= 1 ) THEN
230
231	!
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.
235	p_loc = p
236
237	ENDIF
238
239	!
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!
244	!
245	!-- Left/right
246	IF ( conserve_volume_flow .AND. ( bc_radiation_l .OR. bc_radiation_r ) ) THEN
247
248	volume_flow(1) = 0.0_wp
249	volume_flow_l(1) = 0.0_wp
250
251	IF ( bc_radiation_l ) THEN
252	i = 0
253	ELSEIF ( bc_radiation_r ) THEN
254	i = nx+1
255	ENDIF
256
257	DO j = nys, nyn
258	!
259	!-- Sum up the volume flow through the south/north boundary
260	DO k = nzb+1, nzt
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 ) )
263	ENDDO
264	ENDDO
265
266	#if defined( __parallel )
267	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
268	CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 1, MPI_REAL, MPI_SUM, comm1dy, ierr )
269	#else
270	volume_flow = volume_flow_l
271	#endif
272	volume_flow_offset(1) = ( volume_flow_initial(1) - volume_flow(1) ) / volume_flow_area(1)
273
274	DO j = nysg, nyng
275	DO k = nzb+1, nzt
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 ) )
278	ENDDO
279	ENDDO
280
281	ENDIF
282
283	!
284	!-- South/north
285	IF ( conserve_volume_flow .AND. ( bc_radiation_n .OR. bc_radiation_s ) ) THEN
286
287	volume_flow(2) = 0.0_wp
288	volume_flow_l(2) = 0.0_wp
289
290	IF ( bc_radiation_s ) THEN
291	j = 0
292	ELSEIF ( bc_radiation_n ) THEN
293	j = ny+1
294	ENDIF
295
296	DO i = nxl, nxr
297	!
298	!-- Sum up the volume flow through the south/north boundary
299	DO k = nzb+1, nzt
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 ) )
302	ENDDO
303	ENDDO
304
305	#if defined( __parallel )
306	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
307	CALL MPI_ALLREDUCE( volume_flow_l(2), volume_flow(2), 1, MPI_REAL, MPI_SUM, comm1dx, ierr )
308	#else
309	volume_flow = volume_flow_l
310	#endif
311	volume_flow_offset(2) = ( volume_flow_initial(2) - volume_flow(2) ) / volume_flow_area(2)
312
313	DO i = nxlg, nxrg
314	DO k = nzb+1, nzt
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 ) )
317	ENDDO
318	ENDDO
319
320	ENDIF
321
322	!
323	!-- Remove mean vertical velocity in case that Neumann conditions are used both at bottom and top
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.
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.
333	!-- The removal cannot be done before the first initial time step because ngp_2dh_outer is not yet
334	!-- known then.
335	IF ( ibc_p_b == 1 .AND. ibc_p_t == 1 .AND. .NOT. nesting_offline &
336	.AND. .NOT. ( child_domain .AND. nesting_mode /= 'vertical' ) &
337	.AND. intermediate_timestep_count /= 0 ) THEN
338	w_l = 0.0_wp; w_l_l = 0.0_wp
339	DO i = nxl, nxr
340	DO j = nys, nyn
341	DO k = nzb+1, nzt
342	w_l_l(k) = w_l_l(k) + w(k,j,i) &
343	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
344	ENDDO
345	ENDDO
346	ENDDO
347	#if defined( __parallel )
348	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
349	CALL MPI_ALLREDUCE( w_l_l(1), w_l(1), nzt, MPI_REAL, MPI_SUM, comm2d, ierr )
350	#else
351	w_l = w_l_l
352	#endif
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
358	DO k = nzb+1, nzt
359	w(k,j,i) = w(k,j,i) - w_l(k) &
360	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
361	ENDDO
362	ENDDO
363	ENDDO
364	ENDIF
365
366	!
367	!-- Compute the divergence of the provisional velocity field.
368	CALL cpu_log( log_point_s(1), 'divergence', 'start' )
369
370	IF ( psolver(1:9) == 'multigrid' ) THEN
371	!$OMP PARALLEL DO PRIVATE (i,j,k) SCHEDULE( STATIC )
372	DO i = nxl-1, nxr+1
373	DO j = nys-1, nyn+1
374	DO k = nzb, nzt+1
375	d(k,j,i) = 0.0_wp
376	ENDDO
377	ENDDO
378	ENDDO
379	ELSE
380	!$OMP PARALLEL DO PRIVATE (i,j,k) SCHEDULE( STATIC )
381	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
382	!$ACC PRESENT(d)
383	DO i = nxl, nxr
384	DO j = nys, nyn
385	DO k = nzb+1, nzt
386	d(k,j,i) = 0.0_wp
387	ENDDO
388	ENDDO
389	ENDDO
390	ENDIF
391
392	localsum = 0.0_wp
393	threadsum = 0.0_wp
394
395	!$OMP PARALLEL PRIVATE (i,j,k)
396	!$OMP DO SCHEDULE( STATIC )
397	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
398	!$ACC PRESENT(u, v, w, rho_air, rho_air_zw, ddzw, wall_flags_total_0) &
399	!$ACC PRESENT(d)
400	DO i = nxl, nxr
401	DO j = nys, nyn
402	DO k = 1, nzt
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 ) )
408	ENDDO
409	ENDDO
410	ENDDO
411	!$OMP END PARALLEL
412
413	!
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. &
417	intermediate_timestep_count == 0 ) THEN
418	!$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
419	!$OMP DO SCHEDULE( STATIC )
420	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
421	!$ACC REDUCTION(+:threadsum) COPY(threadsum) &
422	!$ACC PRESENT(d)
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
428	ENDDO
429	ENDDO
430	localsum = localsum + threadsum * dt_3d * weight_pres_l
431	!$OMP END PARALLEL
432	ENDIF
433
434	!
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. &
437	intermediate_timestep_count == 0 ) THEN
438	sums_divold_l(0:statistic_regions) = localsum
439	ENDIF
440
441	CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
442
443	!
444	!-- Compute the pressure perturbation solving the Poisson equation
445	IF ( psolver == 'poisfft' ) THEN
446
447	!
448	!-- Solve Poisson equation via FFT and solution of tridiagonal matrices
449	CALL poisfft( d )
450
451	!
452	!-- Store computed perturbation pressure and set boundary condition in z-direction
453	!$OMP PARALLEL DO PRIVATE (i,j,k) SCHEDULE( STATIC )
454	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
455	!$ACC PRESENT(d, tend)
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:
466	!-- This condition is only required for internal output. The pressure gradient
467	!-- (dp(nzb+1)-dp(nzb))/dz is not used anywhere else.
468	IF ( ibc_p_b == 1 ) THEN
469	!
470	!-- Neumann (dp/dz = 0). Using surface data type, first for non-natural surfaces, then for
471	!-- natural and urban surfaces
472	!-- Upward facing
473	!$OMP PARALLEL DO PRIVATE( m, i, j, k ) SCHEDULE( STATIC )
474	!$ACC PARALLEL LOOP PRIVATE(i, j, k) &
475	!$ACC PRESENT(bc_h, tend)
476	DO m = 1, bc_h(0)%ns
477	i = bc_h(0)%i(m)
478	j = bc_h(0)%j(m)
479	k = bc_h(0)%k(m)
480	tend(k-1,j,i) = tend(k,j,i)
481	ENDDO
482	!
483	!-- Downward facing
484	!$OMP PARALLEL DO PRIVATE( m, i, j, k ) SCHEDULE( STATIC )
485	!$ACC PARALLEL LOOP PRIVATE(i, j, k) &
486	!$ACC PRESENT(bc_h, tend)
487	DO m = 1, bc_h(1)%ns
488	i = bc_h(1)%i(m)
489	j = bc_h(1)%j(m)
490	k = bc_h(1)%k(m)
491	tend(k+1,j,i) = tend(k,j,i)
492	ENDDO
493
494	ELSE
495	!
496	!-- Dirichlet. Using surface data type, first for non-natural surfaces, then for natural and
497	!-- urban surfaces
498	!-- Upward facing
499	!$OMP PARALLEL DO PRIVATE( m, i, j, k ) SCHEDULE( STATIC )
500	DO m = 1, bc_h(0)%ns
501	i = bc_h(0)%i(m)
502	j = bc_h(0)%j(m)
503	k = bc_h(0)%k(m)
504	tend(k-1,j,i) = 0.0_wp
505	ENDDO
506	!
507	!-- Downward facing
508	!$OMP PARALLEL DO PRIVATE( m, i, j, k ) SCHEDULE( STATIC )
509	DO m = 1, bc_h(1)%ns
510	i = bc_h(1)%i(m)
511	j = bc_h(1)%j(m)
512	k = bc_h(1)%k(m)
513	tend(k+1,j,i) = 0.0_wp
514	ENDDO
515
516	ENDIF
517
518	!
519	!-- Top boundary
520	IF ( ibc_p_t == 1 ) THEN
521	!
522	!-- Neumann
523	!$OMP PARALLEL DO PRIVATE (i,j) SCHEDULE( STATIC )
524	DO i = nxlg, nxrg
525	DO j = nysg, nyng
526	tend(nzt+1,j,i) = tend(nzt,j,i)
527	ENDDO
528	ENDDO
529
530	ELSE
531	!
532	!-- Dirichlet
533	!$OMP PARALLEL DO PRIVATE (i,j) SCHEDULE( STATIC )
534	!$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
535	!$ACC PRESENT(tend)
536	DO i = nxlg, nxrg
537	DO j = nysg, nyng
538	tend(nzt+1,j,i) = 0.0_wp
539	ENDDO
540	ENDDO
541
542	ENDIF
543
544	!
545	!-- Exchange boundaries for p
546	CALL exchange_horiz( tend, nbgp )
547
548	ELSEIF ( psolver == 'sor' ) THEN
549
550	!
551	!-- Solve Poisson equation for perturbation pressure using SOR-Red/Black scheme
552	CALL sor( d, ddzu_pres, ddzw, p_loc )
553	tend = p_loc
554
555	ELSEIF ( psolver(1:9) == 'multigrid' ) THEN
556
557	!
558	!-- Solve Poisson equation for perturbation pressure using Multigrid scheme, array tend is used
559	!-- to store the residuals.
560
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.
563	IF ( gathered_size > subdomain_size ) THEN
564	DEALLOCATE( tend )
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) )
568	ENDIF
569
570	IF ( psolver == 'multigrid' ) THEN
571	CALL poismg( tend )
572	ELSE
573	CALL poismg_noopt( tend )
574	ENDIF
575
576	IF ( gathered_size > subdomain_size ) THEN
577	DEALLOCATE( tend )
578	ALLOCATE( tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
579	ENDIF
580
581	!
582	!-- Restore perturbation pressure on tend because this array is used further below to correct the
583	!-- velocity fields
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
591
592	ENDIF
593
594	!
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)
597	IF ( intermediate_timestep_count <= 1 ) THEN
598	!$OMP PARALLEL DO PRIVATE (i,j,k) SCHEDULE( STATIC )
599	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
600	!$ACC PRESENT(p, tend)
601	DO i = nxl-1, nxr+1
602	DO j = nys-1, nyn+1
603	DO k = nzb, nzt+1
604	p(k,j,i) = tend(k,j,i) * weight_substep_l
605	ENDDO
606	ENDDO
607	ENDDO
608
609	ELSEIF ( intermediate_timestep_count > 1 ) THEN
610	!$OMP PARALLEL DO PRIVATE (i,j,k) SCHEDULE( STATIC )
611	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
612	!$ACC PRESENT(p, tend)
613	DO i = nxl-1, nxr+1
614	DO j = nys-1, nyn+1
615	DO k = nzb, nzt+1
616	p(k,j,i) = p(k,j,i) + tend(k,j,i) * weight_substep_l
617	ENDDO
618	ENDDO
619	ENDDO
620
621	ENDIF
622
623	!
624	!-- SOR-method needs ghost layers for the next timestep
625	IF ( psolver == 'sor' ) CALL exchange_horiz( p, nbgp )
626
627	!
628	!-- Correction of the provisional velocities with the current perturbation pressure just computed
629	IF ( conserve_volume_flow .AND. ( bc_lr_cyc .OR. bc_ns_cyc ) ) THEN
630	volume_flow_l(1) = 0.0_wp
631	volume_flow_l(2) = 0.0_wp
632	ENDIF
633	!
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.
639	!$OMP PARALLEL DO PRIVATE (i,j,k) SCHEDULE( STATIC )
640	!$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k) &
641	!$ACC PRESENT(u, v, w, tend, ddzu, wall_flags_total_0)
642	DO i = nxl, nxr
643	DO j = nys, nyn
644
645	DO k = nzb+1, nzt
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 ) )
649	ENDDO
650
651	DO k = nzb+1, nzt
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 ) )
654	ENDDO
655
656	DO k = nzb+1, nzt
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
660
661	ENDDO
662	ENDDO
663
664	!
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
670	w(nzt+1,:,:) = w(nzt,:,:)
671	ENDIF
672
673	!
674	!-- Sum up the volume flow through the right and north boundary
675	IF ( conserve_volume_flow .AND. bc_lr_cyc .AND. bc_ns_cyc .AND. nxr == nx ) THEN
676
677	!$OMP PARALLEL DO PRIVATE (j,k) REDUCTION (+:volume_flow_l(1))
678	DO j = nys, nyn
679	DO k = nzb+1, nzt
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 ) )
682	ENDDO
683	ENDDO
684
685	ENDIF
686
687	IF ( conserve_volume_flow .AND. bc_ns_cyc .AND. bc_lr_cyc .AND. nyn == ny ) THEN
688
689	!$OMP PARALLEL DO PRIVATE (i,k) REDUCTION (+:volume_flow_l(2))
690	DO i = nxl, nxr
691	DO k = nzb+1, nzt
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 ) )
694	ENDDO
695	ENDDO
696
697	ENDIF
698
699	!
700	!-- Conserve the volume flow
701	IF ( conserve_volume_flow .AND. ( bc_lr_cyc .AND. bc_ns_cyc ) ) THEN
702
703	#if defined( __parallel )
704	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
705	CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 2, MPI_REAL, MPI_SUM, comm2d, ierr )
706	#else
707	volume_flow = volume_flow_l
708	#endif
709
710	volume_flow_offset(1:2) = ( volume_flow_initial(1:2) - volume_flow(1:2) ) &
711	/ volume_flow_area(1:2)
712
713	!$OMP PARALLEL DO PRIVATE (i,j,k) SCHEDULE( STATIC )
714	DO i = nxl, nxr
715	DO j = nys, nyn
716	DO k = nzb+1, nzt
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	ENDDO
720	DO k = nzb+1, nzt
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 ) )
723	ENDDO
724	ENDDO
725	ENDDO
726
727	ENDIF
728
729	!
730	!-- Exchange of boundaries for the velocities
731	CALL exchange_horiz( u, nbgp )
732	CALL exchange_horiz( v, nbgp )
733	CALL exchange_horiz( w, nbgp )
734
735	!
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. &
740	intermediate_timestep_count == 0 ) THEN
741	CALL cpu_log( log_point_s(1), 'divergence', 'start' )
742	sums_divnew_l = 0.0_wp
743
744	!
745	!-- d must be reset to zero because it can contain nonzero values below the topography
746	IF ( topography /= 'flat' ) d = 0.0_wp
747
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 )
753	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
754	!$ACC PRESENT(u, v, w, rho_air, rho_air_zw, ddzw, wall_flags_total_0) &
755	!$ACC PRESENT(d)
756	DO i = nxl, nxr
757	DO j = nys, nyn
758	DO k = nzb+1, nzt
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 + &
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 ) )
764	ENDDO
765	ENDDO
766	ENDDO
767	!
768	!-- Compute possible PE-sum of divergences for flow_statistics
769	!$OMP DO SCHEDULE( STATIC )
770	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
771	!$ACC REDUCTION(+:threadsum) COPY(threadsum) &
772	!$ACC PRESENT(d)
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
778	ENDDO
779	ENDDO
780
781	localsum = localsum + threadsum
782	!$OMP END PARALLEL
783
784	!
785	!-- For completeness, set the divergence sum of all statistic regions to those of the total
786	!-- domain
787	sums_divnew_l(0:statistic_regions) = localsum
788
789	CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
790
791	ENDIF
792
793	CALL cpu_log( log_point(8), 'pres', 'stop' )
794
795	END SUBROUTINE pres

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |