Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |