Home

Context Navigation

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

Last change on this file since 3760 was 3655, checked in by knoop, 6 years ago
Bugfix: made "unit" and "found" intend INOUT in module interface subroutines + automatic copyright update
Property svn:keywords set to `Id`
File size: 33.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 3655 2019-01-07 16:51:22Z moh.hefny $
27	! OpenACC port for SPEC
28	!
29	! 3347 2018-10-15 14:21:08Z suehring
30	! Bugfixes in offline nesting.
31	! Add comment.
32	!
33	! 3341 2018-10-15 10:31:27Z suehring
34	! Rename variables for boundary flags and nesting
35	!
36	! 3182 2018-07-27 13:36:03Z suehring
37	! Dollar sign added before Id
38	!
39	! 2696 2017-12-14 17:12:51Z kanani
40	! To avoid jumps while plotting w-profiles w level nzt+1 is set to w level nzt
41	! after velocity modifications through the pressure solver were carried out
42	!
43	! 2696 2017-12-14 17:12:51Z kanani
44	! Corrected "Former revisions" section
45	!
46	! 2696 2017-12-14 17:12:51Z kanani
47	! Change in file header (GPL part)
48	! poismg_noopt modularized (MS)
49	!
50	! 2298 2017-06-29 09:28:18Z raasch
51	! comment changed + some formatting
52	!
53	! 2233 2017-05-30 18:08:54Z suehring
54	!
55	! 2232 2017-05-30 17:47:52Z suehring
56	! Adjustments to new topography and surface concept
57	!
58	! 2118 2017-01-17 16:38:49Z raasch
59	! OpenACC directives and related code removed
60	!
61	! 2073 2016-11-30 14:34:05Z raasch
62	! openmp bugfix for calculation of new divergence
63	!
64	! 2037 2016-10-26 11:15:40Z knoop
65	! Anelastic approximation implemented
66	!
67	! 2000 2016-08-20 18:09:15Z knoop
68	! Forced header and separation lines into 80 columns
69	!
70	! 1932 2016-06-10 12:09:21Z suehring
71	! Initial version of purely vertical nesting introduced.
72	!
73	! 1931 2016-06-10 12:06:59Z suehring
74	! Rename multigrid into multigrid_noopt and multigrid_fast into multigrid
75	!
76	! 1929 2016-06-09 16:25:25Z suehring
77	! Bugfix: weight_substep for initial call, replace by local variable
78	!
79	! 1918 2016-05-27 14:35:57Z raasch
80	! sum of divergence is also calculated when pres is called before the initial
81	! first time step,
82	! stearing is modified by using intermediate_timestep_count = 0 in order to
83	! determine if pres is called before the first initial timestep,
84	! bugfix: in case of Neumann conditions for pressure both at bottom and top,
85	! mean vertical velocity is also removed for the first time step
86	! bugfix for calculating divergences
87	!
88	! 1908 2016-05-25 17:22:32Z suehring
89	! New divergence for output into RUN_CONTROL file is calculated only at last
90	! Runge-Kutta step
91	!
92	! 1845 2016-04-08 08:29:13Z raasch
93	! nzb_2d replace by nzb_u\|v_inner
94	!
95	! 1799 2016-04-05 08:35:55Z gronemeier
96	! Bugfix: excluded third dimension from horizontal volume flow calculation
97	!
98	! 1762 2016-02-25 12:31:13Z hellstea
99	! Introduction of nested domain feature
100	!
101	! 1682 2015-10-07 23:56:08Z knoop
102	! Code annotations made doxygen readable
103	!
104	! 1575 2015-03-27 09:56:27Z raasch
105	! poismg_fast + respective module added, adjustments for psolver-queries
106	!
107	! 1342 2014-03-26 17:04:47Z kanani
108	! REAL constants defined as wp-kind
109	!
110	! 1320 2014-03-20 08:40:49Z raasch
111	! ONLY-attribute added to USE-statements,
112	! kind-parameters added to all INTEGER and REAL declaration statements,
113	! kinds are defined in new module kinds,
114	! old module precision_kind is removed,
115	! revision history before 2012 removed,
116	! comment fields (!:) to be used for variable explanations added to
117	! all variable declaration statements
118	!
119	! 1318 2014-03-17 13:35:16Z raasch
120	! module interfaces removed
121	!
122	! 1306 2014-03-13 14:30:59Z raasch
123	! second argument removed from routine poisfft
124	!
125	! 1257 2013-11-08 15:18:40Z raasch
126	! openacc loop and loop vector clauses removed, independent clauses added,
127	! end parallel replaced by end parallel loop
128	!
129	! 1221 2013-09-10 08:59:13Z raasch
130	! openACC porting of reduction operations, loops for calculating d are
131	! using the rflags_s_inner multiply flag instead of the nzb_s_inner loop index
132	!
133	! 1212 2013-08-15 08:46:27Z raasch
134	! call of poisfft_hybrid removed
135	!
136	! 1117 2013-03-27 11:15:36Z suehring
137	! Bugfix in OpenMP parallelization.
138	!
139	! 1113 2013-03-10 02:48:14Z raasch
140	! GPU-porting of several loops, some loops rearranged
141	!
142	! 1111 2013-03-08 23:54:10Z
143	! openACC statements added,
144	! ibc_p_b = 2 removed
145	!
146	! 1092 2013-02-02 11:24:22Z raasch
147	! unused variables removed
148	!
149	! 1036 2012-10-22 13:43:42Z raasch
150	! code put under GPL (PALM 3.9)
151	!
152	! 1003 2012-09-14 14:35:53Z raasch
153	! adjustment of array tend for cases with unequal subdomain sizes removed
154	!
155	! Revision 1.1 1997/07/24 11:24:44 raasch
156	! Initial revision
157	!
158	!
159	! Description:
160	! ------------
161	!> Compute the divergence of the provisional velocity field. Solve the Poisson
162	!> equation for the perturbation pressure. Compute the final velocities using
163	!> this perturbation pressure. Compute the remaining divergence.
164	!------------------------------------------------------------------------------!
165	SUBROUTINE pres
166
167
168	USE arrays_3d, &
169	ONLY: d, ddzu, ddzu_pres, ddzw, dzw, p, p_loc, rho_air, rho_air_zw, &
170	tend, u, v, w
171
172	USE control_parameters, &
173	ONLY: bc_lr_cyc, bc_ns_cyc, bc_radiation_l, bc_radiation_n, &
174	bc_radiation_r, bc_radiation_s, child_domain, &
175	conserve_volume_flow, coupling_mode, &
176	dt_3d, gathered_size, ibc_p_b, ibc_p_t, &
177	intermediate_timestep_count, intermediate_timestep_count_max, &
178	mg_switch_to_pe0_level, nesting_offline, &
179	psolver, subdomain_size, &
180	topography, volume_flow, volume_flow_area, volume_flow_initial
181
182	USE cpulog, &
183	ONLY: cpu_log, log_point, log_point_s
184
185	USE grid_variables, &
186	ONLY: ddx, ddy
187
188	USE indices, &
189	ONLY: nbgp, ngp_2dh_outer, nx, nxl, nxlg, nxl_mg, nxr, nxrg, nxr_mg, &
190	ny, nys, nysg, nys_mg, nyn, nyng, nyn_mg, nzb, nzt, nzt_mg, &
191	wall_flags_0
192
193	USE kinds
194
195	USE pegrid
196
197	USE pmc_interface, &
198	ONLY: nesting_mode
199
200	USE poisfft_mod, &
201	ONLY: poisfft
202
203	USE poismg_mod
204
205	USE poismg_noopt_mod
206
207	USE statistics, &
208	ONLY: statistic_regions, sums_divnew_l, sums_divold_l, weight_pres, &
209	weight_substep
210
211	USE surface_mod, &
212	ONLY : bc_h
213
214	IMPLICIT NONE
215
216	INTEGER(iwp) :: i !<
217	INTEGER(iwp) :: j !<
218	INTEGER(iwp) :: k !<
219	INTEGER(iwp) :: m !<
220
221	REAL(wp) :: ddt_3d !<
222	REAL(wp) :: d_weight_pres !<
223	REAL(wp) :: localsum !<
224	REAL(wp) :: threadsum !<
225	REAL(wp) :: weight_pres_l !<
226	REAL(wp) :: weight_substep_l !<
227
228	REAL(wp), DIMENSION(1:3) :: volume_flow_l !<
229	REAL(wp), DIMENSION(1:3) :: volume_flow_offset !<
230	REAL(wp), DIMENSION(1:nzt) :: w_l !<
231	REAL(wp), DIMENSION(1:nzt) :: w_l_l !<
232
233	LOGICAL :: child_domain_nvn !<
234
235
236	CALL cpu_log( log_point(8), 'pres', 'start' )
237
238	!$ACC DATA COPYOUT(d, tend) &
239	!$ACC COPY(u, v, w, p) &
240	!$ACC COPYIN(rho_air, rho_air_zw, ddzu, ddzw, wall_flags_0)
241
242	!$ACC DATA &
243	!$ACC COPYIN(bc_h(0:1)) &
244	!$ACC COPYIN(bc_h(0)%i(1:bc_h(0)%ns)) &
245	!$ACC COPYIN(bc_h(0)%j(1:bc_h(0)%ns)) &
246	!$ACC COPYIN(bc_h(0)%k(1:bc_h(0)%ns)) &
247	!$ACC COPYIN(bc_h(1)%i(1:bc_h(1)%ns)) &
248	!$ACC COPYIN(bc_h(1)%j(1:bc_h(1)%ns)) &
249	!$ACC COPYIN(bc_h(1)%k(1:bc_h(1)%ns))
250
251
252	!
253	!-- Calculate quantities to be used locally
254	ddt_3d = 1.0_wp / dt_3d
255	IF ( intermediate_timestep_count == 0 ) THEN
256	!
257	!-- If pres is called before initial time step
258	weight_pres_l = 1.0_wp
259	d_weight_pres = 1.0_wp
260	weight_substep_l = 1.0_wp
261	ELSE
262	weight_pres_l = weight_pres(intermediate_timestep_count)
263	d_weight_pres = 1.0_wp / weight_pres(intermediate_timestep_count)
264	weight_substep_l = weight_substep(intermediate_timestep_count)
265	ENDIF
266
267	!
268	!-- Multigrid method expects array d to have one ghost layer.
269	!--
270	IF ( psolver(1:9) == 'multigrid' ) THEN
271
272	DEALLOCATE( d )
273	ALLOCATE( d(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
274
275	!
276	!-- Since p is later used to hold the weighted average of the substeps, it
277	!-- cannot be used in the iterative solver. Therefore, its initial value is
278	!-- stored on p_loc, which is then iteratively advanced in every substep.
279	IF ( intermediate_timestep_count <= 1 ) THEN
280	DO i = nxl-1, nxr+1
281	DO j = nys-1, nyn+1
282	DO k = nzb, nzt+1
283	p_loc(k,j,i) = p(k,j,i)
284	ENDDO
285	ENDDO
286	ENDDO
287	ENDIF
288
289	ELSEIF ( psolver == 'sor' .AND. intermediate_timestep_count <= 1 ) THEN
290
291	!
292	!-- Since p is later used to hold the weighted average of the substeps, it
293	!-- cannot be used in the iterative solver. Therefore, its initial value is
294	!-- stored on p_loc, which is then iteratively advanced in every substep.
295	p_loc = p
296
297	ENDIF
298
299	!
300	!-- Conserve the volume flow at the outflow in case of non-cyclic lateral
301	!-- boundary conditions
302	!-- WARNING: so far, this conservation does not work at the left/south
303	!-- boundary if the topography at the inflow differs from that at the
304	!-- outflow! For this case, volume_flow_area needs adjustment!
305	!
306	!-- Left/right
307	IF ( conserve_volume_flow .AND. ( bc_radiation_l .OR. &
308	bc_radiation_r ) ) THEN
309
310	volume_flow(1) = 0.0_wp
311	volume_flow_l(1) = 0.0_wp
312
313	IF ( bc_radiation_l ) THEN
314	i = 0
315	ELSEIF ( bc_radiation_r ) THEN
316	i = nx+1
317	ENDIF
318
319	DO j = nys, nyn
320	!
321	!-- Sum up the volume flow through the south/north boundary
322	DO k = nzb+1, nzt
323	volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k) &
324	* MERGE( 1.0_wp, 0.0_wp, &
325	BTEST( wall_flags_0(k,j,i), 1 ) &
326	)
327	ENDDO
328	ENDDO
329
330	#if defined( __parallel )
331	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
332	CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 1, MPI_REAL, &
333	MPI_SUM, comm1dy, ierr )
334	#else
335	volume_flow = volume_flow_l
336	#endif
337	volume_flow_offset(1) = ( volume_flow_initial(1) - volume_flow(1) ) &
338	/ volume_flow_area(1)
339
340	DO j = nysg, nyng
341	DO k = nzb+1, nzt
342	u(k,j,i) = u(k,j,i) + volume_flow_offset(1) &
343	* MERGE( 1.0_wp, 0.0_wp, &
344	BTEST( wall_flags_0(k,j,i), 1 ) &
345	)
346	ENDDO
347	ENDDO
348
349	ENDIF
350
351	!
352	!-- South/north
353	IF ( conserve_volume_flow .AND. ( bc_radiation_n .OR. bc_radiation_s ) ) THEN
354
355	volume_flow(2) = 0.0_wp
356	volume_flow_l(2) = 0.0_wp
357
358	IF ( bc_radiation_s ) THEN
359	j = 0
360	ELSEIF ( bc_radiation_n ) THEN
361	j = ny+1
362	ENDIF
363
364	DO i = nxl, nxr
365	!
366	!-- Sum up the volume flow through the south/north boundary
367	DO k = nzb+1, nzt
368	volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k) &
369	* MERGE( 1.0_wp, 0.0_wp, &
370	BTEST( wall_flags_0(k,j,i), 2 ) &
371	)
372	ENDDO
373	ENDDO
374
375	#if defined( __parallel )
376	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
377	CALL MPI_ALLREDUCE( volume_flow_l(2), volume_flow(2), 1, MPI_REAL, &
378	MPI_SUM, comm1dx, ierr )
379	#else
380	volume_flow = volume_flow_l
381	#endif
382	volume_flow_offset(2) = ( volume_flow_initial(2) - volume_flow(2) ) &
383	/ volume_flow_area(2)
384
385	DO i = nxlg, nxrg
386	DO k = nzb+1, nzt
387	v(k,j,i) = v(k,j,i) + volume_flow_offset(2) &
388	* MERGE( 1.0_wp, 0.0_wp, &
389	BTEST( wall_flags_0(k,j,i), 2 ) &
390	)
391	ENDDO
392	ENDDO
393
394	ENDIF
395
396	!
397	!-- Remove mean vertical velocity in case that Neumann conditions are
398	!-- used both at bottom and top boundary, and if not a nested domain in a
399	!-- normal nesting run. In case of vertical nesting, this must be done.
400	!-- Therefore an auxiliary logical variable child_domain_nvn is used here, and
401	!-- nvn stands for non-vertical nesting.
402	!-- This cannot be done before the first initial time step because ngp_2dh_outer
403	!-- is not yet known then.
404	child_domain_nvn = child_domain
405	IF ( child_domain .AND. nesting_mode == 'vertical' ) THEN
406	child_domain_nvn = .FALSE.
407	ENDIF
408
409	IF ( ibc_p_b == 1 .AND. ibc_p_t == 1 .AND. .NOT. nesting_offline .AND. &
410	.NOT. child_domain_nvn .AND. intermediate_timestep_count /= 0 ) &
411	THEN
412	w_l = 0.0_wp; w_l_l = 0.0_wp
413	DO i = nxl, nxr
414	DO j = nys, nyn
415	DO k = nzb+1, nzt
416	w_l_l(k) = w_l_l(k) + w(k,j,i) &
417	* MERGE( 1.0_wp, 0.0_wp, &
418	BTEST( wall_flags_0(k,j,i), 3 ) &
419	)
420	ENDDO
421	ENDDO
422	ENDDO
423	#if defined( __parallel )
424	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
425	CALL MPI_ALLREDUCE( w_l_l(1), w_l(1), nzt, MPI_REAL, MPI_SUM, &
426	comm2d, ierr )
427	#else
428	w_l = w_l_l
429	#endif
430	DO k = 1, nzt
431	w_l(k) = w_l(k) / ngp_2dh_outer(k,0)
432	ENDDO
433	DO i = nxlg, nxrg
434	DO j = nysg, nyng
435	DO k = nzb+1, nzt
436	w(k,j,i) = w(k,j,i) - w_l(k) &
437	* MERGE( 1.0_wp, 0.0_wp, &
438	BTEST( wall_flags_0(k,j,i), 3 ) &
439	)
440	ENDDO
441	ENDDO
442	ENDDO
443	ENDIF
444
445	!
446	!-- Compute the divergence of the provisional velocity field.
447	CALL cpu_log( log_point_s(1), 'divergence', 'start' )
448
449	IF ( psolver(1:9) == 'multigrid' ) THEN
450	!$OMP PARALLEL DO SCHEDULE( STATIC ) PRIVATE (i,j,k)
451	DO i = nxl-1, nxr+1
452	DO j = nys-1, nyn+1
453	DO k = nzb, nzt+1
454	d(k,j,i) = 0.0_wp
455	ENDDO
456	ENDDO
457	ENDDO
458	ELSE
459	!$OMP PARALLEL DO SCHEDULE( STATIC ) PRIVATE (i,j,k)
460	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
461	!$ACC PRESENT(d)
462	DO i = nxl, nxr
463	DO j = nys, nyn
464	DO k = nzb+1, nzt
465	d(k,j,i) = 0.0_wp
466	ENDDO
467	ENDDO
468	ENDDO
469	ENDIF
470
471	localsum = 0.0_wp
472	threadsum = 0.0_wp
473
474	#if defined( __ibm )
475	!$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
476	!$OMP DO SCHEDULE( STATIC )
477	DO i = nxl, nxr
478	DO j = nys, nyn
479	DO k = nzb+1, nzt
480	d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx + &
481	( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy + &
482	( w(k,j,i) * rho_air_zw(k) - &
483	w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k) &
484	) * ddt_3d * d_weight_pres &
485	* MERGE( 1.0_wp, 0.0_wp, &
486	BTEST( wall_flags_0(k,j,i), 0 ) &
487	)
488	ENDDO
489	!
490	!-- Compute possible PE-sum of divergences for flow_statistics
491	DO k = nzb+1, nzt
492	threadsum = threadsum + ABS( d(k,j,i) ) &
493	* MERGE( 1.0_wp, 0.0_wp, &
494	BTEST( wall_flags_0(k,j,i), 0 ) &
495	)
496	ENDDO
497
498	ENDDO
499	ENDDO
500
501	IF ( intermediate_timestep_count == intermediate_timestep_count_max .OR. &
502	intermediate_timestep_count == 0 ) THEN
503	localsum = localsum + threadsum * dt_3d * weight_pres_l
504	ENDIF
505	!$OMP END PARALLEL
506	#else
507
508	!$OMP PARALLEL PRIVATE (i,j,k)
509	!$OMP DO SCHEDULE( STATIC )
510	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
511	!$ACC PRESENT(u, v, w, rho_air, rho_air_zw, ddzw, wall_flags_0) &
512	!$ACC PRESENT(d)
513	DO i = nxl, nxr
514	DO j = nys, nyn
515	DO k = 1, nzt
516	d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx + &
517	( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy + &
518	( w(k,j,i) * rho_air_zw(k) - &
519	w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k) &
520	) * ddt_3d * d_weight_pres &
521	* MERGE( 1.0_wp, 0.0_wp, &
522	BTEST( wall_flags_0(k,j,i), 0 ) &
523	)
524	ENDDO
525	ENDDO
526	ENDDO
527	!$OMP END PARALLEL
528
529	!
530	!-- Compute possible PE-sum of divergences for flow_statistics. Carry out
531	!-- computation only at last Runge-Kutta substep.
532	IF ( intermediate_timestep_count == intermediate_timestep_count_max .OR. &
533	intermediate_timestep_count == 0 ) THEN
534	!$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
535	!$OMP DO SCHEDULE( STATIC )
536	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
537	!$ACC REDUCTION(+:threadsum) COPY(threadsum) &
538	!$ACC PRESENT(d)
539	DO i = nxl, nxr
540	DO j = nys, nyn
541	DO k = nzb+1, nzt
542	threadsum = threadsum + ABS( d(k,j,i) )
543	ENDDO
544	ENDDO
545	ENDDO
546	localsum = localsum + threadsum * dt_3d * weight_pres_l
547	!$OMP END PARALLEL
548	ENDIF
549	#endif
550
551	!
552	!-- For completeness, set the divergence sum of all statistic regions to those
553	!-- of the total domain
554	IF ( intermediate_timestep_count == intermediate_timestep_count_max .OR. &
555	intermediate_timestep_count == 0 ) THEN
556	sums_divold_l(0:statistic_regions) = localsum
557	ENDIF
558
559	CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
560
561	!
562	!-- Compute the pressure perturbation solving the Poisson equation
563	IF ( psolver == 'poisfft' ) THEN
564
565	!
566	!-- Solve Poisson equation via FFT and solution of tridiagonal matrices
567	CALL poisfft( d )
568
569	!
570	!-- Store computed perturbation pressure and set boundary condition in
571	!-- z-direction
572	!$OMP PARALLEL DO PRIVATE (i,j,k)
573	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
574	!$ACC PRESENT(d, tend)
575	DO i = nxl, nxr
576	DO j = nys, nyn
577	DO k = nzb+1, nzt
578	tend(k,j,i) = d(k,j,i)
579	ENDDO
580	ENDDO
581	ENDDO
582
583	!
584	!-- Bottom boundary:
585	!-- This condition is only required for internal output. The pressure
586	!-- gradient (dp(nzb+1)-dp(nzb))/dz is not used anywhere else.
587	IF ( ibc_p_b == 1 ) THEN
588	!
589	!-- Neumann (dp/dz = 0). Using surfae data type, first for non-natural
590	!-- surfaces, then for natural and urban surfaces
591	!-- Upward facing
592	!$OMP PARALLEL DO PRIVATE( i, j, k )
593	!$ACC PARALLEL LOOP PRIVATE(i, j, k) &
594	!$ACC PRESENT(bc_h, tend)
595	DO m = 1, bc_h(0)%ns
596	i = bc_h(0)%i(m)
597	j = bc_h(0)%j(m)
598	k = bc_h(0)%k(m)
599	tend(k-1,j,i) = tend(k,j,i)
600	ENDDO
601	!
602	!-- Downward facing
603	!$OMP PARALLEL DO PRIVATE( i, j, k )
604	!$ACC PARALLEL LOOP PRIVATE(i, j, k) &
605	!$ACC PRESENT(bc_h, tend)
606	DO m = 1, bc_h(1)%ns
607	i = bc_h(1)%i(m)
608	j = bc_h(1)%j(m)
609	k = bc_h(1)%k(m)
610	tend(k+1,j,i) = tend(k,j,i)
611	ENDDO
612
613	ELSE
614	!
615	!-- Dirichlet. Using surface data type, first for non-natural
616	!-- surfaces, then for natural and urban surfaces
617	!-- Upward facing
618	!$OMP PARALLEL DO PRIVATE( i, j, k )
619	DO m = 1, bc_h(0)%ns
620	i = bc_h(0)%i(m)
621	j = bc_h(0)%j(m)
622	k = bc_h(0)%k(m)
623	tend(k-1,j,i) = 0.0_wp
624	ENDDO
625	!
626	!-- Downward facing
627	!$OMP PARALLEL DO PRIVATE( i, j, k )
628	DO m = 1, bc_h(1)%ns
629	i = bc_h(1)%i(m)
630	j = bc_h(1)%j(m)
631	k = bc_h(1)%k(m)
632	tend(k+1,j,i) = 0.0_wp
633	ENDDO
634
635	ENDIF
636
637	!
638	!-- Top boundary
639	IF ( ibc_p_t == 1 ) THEN
640	!
641	!-- Neumann
642	!$OMP PARALLEL DO PRIVATE (i,j,k)
643	DO i = nxlg, nxrg
644	DO j = nysg, nyng
645	tend(nzt+1,j,i) = tend(nzt,j,i)
646	ENDDO
647	ENDDO
648
649	ELSE
650	!
651	!-- Dirichlet
652	!$OMP PARALLEL DO PRIVATE (i,j,k)
653	!$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
654	!$ACC PRESENT(tend)
655	DO i = nxlg, nxrg
656	DO j = nysg, nyng
657	tend(nzt+1,j,i) = 0.0_wp
658	ENDDO
659	ENDDO
660
661	ENDIF
662
663	!
664	!-- Exchange boundaries for p
665	CALL exchange_horiz( tend, nbgp )
666
667	ELSEIF ( psolver == 'sor' ) THEN
668
669	!
670	!-- Solve Poisson equation for perturbation pressure using SOR-Red/Black
671	!-- scheme
672	CALL sor( d, ddzu_pres, ddzw, p_loc )
673	tend = p_loc
674
675	ELSEIF ( psolver(1:9) == 'multigrid' ) THEN
676
677	!
678	!-- Solve Poisson equation for perturbation pressure using Multigrid scheme,
679	!-- array tend is used to store the residuals.
680
681	!-- If the number of grid points of the gathered grid, which is collected
682	!-- on PE0, is larger than the number of grid points of an PE, than array
683	!-- tend will be enlarged.
684	IF ( gathered_size > subdomain_size ) THEN
685	DEALLOCATE( tend )
686	ALLOCATE( tend(nzb:nzt_mg(mg_switch_to_pe0_level)+1,nys_mg( &
687	mg_switch_to_pe0_level)-1:nyn_mg(mg_switch_to_pe0_level)+1,&
688	nxl_mg(mg_switch_to_pe0_level)-1:nxr_mg( &
689	mg_switch_to_pe0_level)+1) )
690	ENDIF
691
692	IF ( psolver == 'multigrid' ) THEN
693	CALL poismg( tend )
694	ELSE
695	CALL poismg_noopt( tend )
696	ENDIF
697
698	IF ( gathered_size > subdomain_size ) THEN
699	DEALLOCATE( tend )
700	ALLOCATE( tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
701	ENDIF
702
703	!
704	!-- Restore perturbation pressure on tend because this array is used
705	!-- further below to correct the velocity fields
706	DO i = nxl-1, nxr+1
707	DO j = nys-1, nyn+1
708	DO k = nzb, nzt+1
709	tend(k,j,i) = p_loc(k,j,i)
710	ENDDO
711	ENDDO
712	ENDDO
713
714	ENDIF
715
716	!
717	!-- Store perturbation pressure on array p, used for pressure data output.
718	!-- Ghost layers are added in the output routines (except sor-method: see below)
719	IF ( intermediate_timestep_count <= 1 ) THEN
720	!$OMP PARALLEL PRIVATE (i,j,k)
721	!$OMP DO
722	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
723	!$ACC PRESENT(p, tend)
724	DO i = nxl-1, nxr+1
725	DO j = nys-1, nyn+1
726	DO k = nzb, nzt+1
727	p(k,j,i) = tend(k,j,i) * &
728	weight_substep_l
729	ENDDO
730	ENDDO
731	ENDDO
732	!$OMP END PARALLEL
733
734	ELSEIF ( intermediate_timestep_count > 1 ) THEN
735	!$OMP PARALLEL PRIVATE (i,j,k)
736	!$OMP DO
737	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
738	!$ACC PRESENT(p, tend)
739	DO i = nxl-1, nxr+1
740	DO j = nys-1, nyn+1
741	DO k = nzb, nzt+1
742	p(k,j,i) = p(k,j,i) + tend(k,j,i) * &
743	weight_substep_l
744	ENDDO
745	ENDDO
746	ENDDO
747	!$OMP END PARALLEL
748
749	ENDIF
750
751	!
752	!-- SOR-method needs ghost layers for the next timestep
753	IF ( psolver == 'sor' ) CALL exchange_horiz( p, nbgp )
754
755	!
756	!-- Correction of the provisional velocities with the current perturbation
757	!-- pressure just computed
758	IF ( conserve_volume_flow .AND. ( bc_lr_cyc .OR. bc_ns_cyc ) ) THEN
759	volume_flow_l(1) = 0.0_wp
760	volume_flow_l(2) = 0.0_wp
761	ENDIF
762	!
763	!-- Add pressure gradients to the velocity components. Note, the loops are
764	!-- running over the entire model domain, even though, in case of non-cyclic
765	!-- boundaries u- and v-component are not prognostic at i=0 and j=0,
766	!-- respectiveley. However, in case of Dirichlet boundary conditions for the
767	!-- velocities, zero-gradient conditions for the pressure are set, so that
768	!-- no modification is imposed at the boundaries.
769	!$OMP PARALLEL PRIVATE (i,j,k)
770	!$OMP DO
771	!$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k) &
772	!$ACC PRESENT(u, v, w, tend, ddzu, wall_flags_0)
773	DO i = nxl, nxr
774	DO j = nys, nyn
775
776	DO k = nzb+1, nzt
777	w(k,j,i) = w(k,j,i) - dt_3d * &
778	( tend(k+1,j,i) - tend(k,j,i) ) * ddzu(k+1) &
779	* weight_pres_l &
780	* MERGE( 1.0_wp, 0.0_wp, &
781	BTEST( wall_flags_0(k,j,i), 3 ) &
782	)
783	ENDDO
784
785	DO k = nzb+1, nzt
786	u(k,j,i) = u(k,j,i) - dt_3d * &
787	( tend(k,j,i) - tend(k,j,i-1) ) * ddx &
788	* weight_pres_l &
789	* MERGE( 1.0_wp, 0.0_wp, &
790	BTEST( wall_flags_0(k,j,i), 1 ) &
791	)
792	ENDDO
793
794	DO k = nzb+1, nzt
795	v(k,j,i) = v(k,j,i) - dt_3d * &
796	( tend(k,j,i) - tend(k,j-1,i) ) * ddy &
797	* weight_pres_l &
798	* MERGE( 1.0_wp, 0.0_wp, &
799	BTEST( wall_flags_0(k,j,i), 2 ) &
800	)
801	ENDDO
802
803	ENDDO
804	ENDDO
805	!$OMP END PARALLEL
806
807	!
808	!-- The vertical velocity is not set to zero at nzt + 1 for nested domains
809	!-- Instead it is set to the values of nzt (see routine vnest_boundary_conds
810	!-- or pmci_interp_tril_t) BEFORE calling the pressure solver. To avoid jumps
811	!-- while plotting profiles w at the top has to be set to the values in the
812	!-- height nzt after above modifications. Hint: w level nzt+1 does not impact
813	!-- results.
814	IF ( child_domain .OR. coupling_mode == 'vnested_fine' ) THEN
815	w(nzt+1,:,:) = w(nzt,:,:)
816	ENDIF
817
818	!
819	!-- Sum up the volume flow through the right and north boundary
820	IF ( conserve_volume_flow .AND. bc_lr_cyc .AND. bc_ns_cyc .AND. &
821	nxr == nx ) THEN
822
823	!$OMP PARALLEL PRIVATE (j,k)
824	!$OMP DO
825	DO j = nys, nyn
826	!$OMP CRITICAL
827	DO k = nzb+1, nzt
828	volume_flow_l(1) = volume_flow_l(1) + u(k,j,nxr) * dzw(k) &
829	* MERGE( 1.0_wp, 0.0_wp, &
830	BTEST( wall_flags_0(k,j,nxr), 1 )&
831	)
832	ENDDO
833	!$OMP END CRITICAL
834	ENDDO
835	!$OMP END PARALLEL
836
837	ENDIF
838
839	IF ( conserve_volume_flow .AND. bc_ns_cyc .AND. bc_lr_cyc .AND. &
840	nyn == ny ) THEN
841
842	!$OMP PARALLEL PRIVATE (i,k)
843	!$OMP DO
844	DO i = nxl, nxr
845	!$OMP CRITICAL
846	DO k = nzb+1, nzt
847	volume_flow_l(2) = volume_flow_l(2) + v(k,nyn,i) * dzw(k) &
848	* MERGE( 1.0_wp, 0.0_wp, &
849	BTEST( wall_flags_0(k,nyn,i), 2 )&
850	)
851	ENDDO
852	!$OMP END CRITICAL
853	ENDDO
854	!$OMP END PARALLEL
855
856	ENDIF
857
858	!
859	!-- Conserve the volume flow
860	IF ( conserve_volume_flow .AND. ( bc_lr_cyc .AND. bc_ns_cyc ) ) THEN
861
862	#if defined( __parallel )
863	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
864	CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 2, MPI_REAL, &
865	MPI_SUM, comm2d, ierr )
866	#else
867	volume_flow = volume_flow_l
868	#endif
869
870	volume_flow_offset(1:2) = ( volume_flow_initial(1:2) - volume_flow(1:2) ) / &
871	volume_flow_area(1:2)
872
873	!$OMP PARALLEL PRIVATE (i,j,k)
874	!$OMP DO
875	DO i = nxl, nxr
876	DO j = nys, nyn
877	DO k = nzb+1, nzt
878	u(k,j,i) = u(k,j,i) + volume_flow_offset(1) &
879	* MERGE( 1.0_wp, 0.0_wp, &
880	BTEST( wall_flags_0(k,j,i), 1 ) &
881	)
882	ENDDO
883	DO k = nzb+1, nzt
884	v(k,j,i) = v(k,j,i) + volume_flow_offset(2) &
885	* MERGE( 1.0_wp, 0.0_wp, &
886	BTEST( wall_flags_0(k,j,i), 2 ) &
887	)
888	ENDDO
889	ENDDO
890	ENDDO
891
892	!$OMP END PARALLEL
893
894	ENDIF
895
896	!
897	!-- Exchange of boundaries for the velocities
898	CALL exchange_horiz( u, nbgp )
899	CALL exchange_horiz( v, nbgp )
900	CALL exchange_horiz( w, nbgp )
901
902	!
903	!-- Compute the divergence of the corrected velocity field,
904	!-- A possible PE-sum is computed in flow_statistics. Carry out computation
905	!-- only at last Runge-Kutta step.
906	IF ( intermediate_timestep_count == intermediate_timestep_count_max .OR. &
907	intermediate_timestep_count == 0 ) THEN
908	CALL cpu_log( log_point_s(1), 'divergence', 'start' )
909	sums_divnew_l = 0.0_wp
910
911	!
912	!-- d must be reset to zero because it can contain nonzero values below the
913	!-- topography
914	IF ( topography /= 'flat' ) d = 0.0_wp
915
916	localsum = 0.0_wp
917	threadsum = 0.0_wp
918
919	!$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
920	#if defined( __ibm )
921	!$OMP DO SCHEDULE( STATIC )
922	DO i = nxl, nxr
923	DO j = nys, nyn
924	DO k = nzb+1, nzt
925	d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx + &
926	( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy + &
927	( w(k,j,i) * rho_air_zw(k) - &
928	w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k) &
929	) * MERGE( 1.0_wp, 0.0_wp, &
930	BTEST( wall_flags_0(k,j,i), 0 ) &
931	)
932	ENDDO
933	DO k = nzb+1, nzt
934	threadsum = threadsum + ABS( d(k,j,i) )
935	ENDDO
936	ENDDO
937	ENDDO
938	#else
939	!$OMP DO SCHEDULE( STATIC )
940	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
941	!$ACC PRESENT(u, v, w, rho_air, rho_air_zw, ddzw, wall_flags_0) &
942	!$ACC PRESENT(d)
943	DO i = nxl, nxr
944	DO j = nys, nyn
945	DO k = nzb+1, nzt
946	d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx + &
947	( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy + &
948	( w(k,j,i) * rho_air_zw(k) - &
949	w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k) &
950	) * MERGE( 1.0_wp, 0.0_wp, &
951	BTEST( wall_flags_0(k,j,i), 0 ) &
952	)
953	ENDDO
954	ENDDO
955	ENDDO
956	!
957	!-- Compute possible PE-sum of divergences for flow_statistics
958	!$OMP DO SCHEDULE( STATIC )
959	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
960	!$ACC REDUCTION(+:threadsum) COPY(threadsum) &
961	!$ACC PRESENT(d)
962	DO i = nxl, nxr
963	DO j = nys, nyn
964	DO k = nzb+1, nzt
965	threadsum = threadsum + ABS( d(k,j,i) )
966	ENDDO
967	ENDDO
968	ENDDO
969	#endif
970
971	localsum = localsum + threadsum
972	!$OMP END PARALLEL
973
974	!
975	!-- For completeness, set the divergence sum of all statistic regions to those
976	!-- of the total domain
977	sums_divnew_l(0:statistic_regions) = localsum
978
979	CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
980
981	ENDIF
982
983	CALL cpu_log( log_point(8), 'pres', 'stop' )
984
985	!$ACC END DATA
986	!$ACC END DATA
987
988	END SUBROUTINE pres

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |