Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |