Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |