Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |