Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |