Home

Context Navigation

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

Last change on this file since 3949 was 3849, checked in by knoop, 6 years ago
Bugfix: added proper OpenACC support to disturb_field
Property svn:keywords set to `Id`
File size: 33.5 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 3849 2019-04-01 16:35:16Z forkel $
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	!
239	!-- Calculate quantities to be used locally
240	ddt_3d = 1.0_wp / dt_3d
241	IF ( intermediate_timestep_count == 0 ) THEN
242	!
243	!-- If pres is called before initial time step
244	weight_pres_l = 1.0_wp
245	d_weight_pres = 1.0_wp
246	weight_substep_l = 1.0_wp
247	ELSE
248	weight_pres_l = weight_pres(intermediate_timestep_count)
249	d_weight_pres = 1.0_wp / weight_pres(intermediate_timestep_count)
250	weight_substep_l = weight_substep(intermediate_timestep_count)
251	ENDIF
252
253	!
254	!-- Multigrid method expects array d to have one ghost layer.
255	!--
256	IF ( psolver(1:9) == 'multigrid' ) THEN
257
258	DEALLOCATE( d )
259	ALLOCATE( d(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
260
261	!
262	!-- Since p is later used to hold the weighted average of the substeps, it
263	!-- cannot be used in the iterative solver. Therefore, its initial value is
264	!-- stored on p_loc, which is then iteratively advanced in every substep.
265	IF ( intermediate_timestep_count <= 1 ) THEN
266	DO i = nxl-1, nxr+1
267	DO j = nys-1, nyn+1
268	DO k = nzb, nzt+1
269	p_loc(k,j,i) = p(k,j,i)
270	ENDDO
271	ENDDO
272	ENDDO
273	ENDIF
274
275	ELSEIF ( psolver == 'sor' .AND. intermediate_timestep_count <= 1 ) THEN
276
277	!
278	!-- Since p is later used to hold the weighted average of the substeps, it
279	!-- cannot be used in the iterative solver. Therefore, its initial value is
280	!-- stored on p_loc, which is then iteratively advanced in every substep.
281	p_loc = p
282
283	ENDIF
284
285	!
286	!-- Conserve the volume flow at the outflow in case of non-cyclic lateral
287	!-- boundary conditions
288	!-- WARNING: so far, this conservation does not work at the left/south
289	!-- boundary if the topography at the inflow differs from that at the
290	!-- outflow! For this case, volume_flow_area needs adjustment!
291	!
292	!-- Left/right
293	IF ( conserve_volume_flow .AND. ( bc_radiation_l .OR. &
294	bc_radiation_r ) ) THEN
295
296	volume_flow(1) = 0.0_wp
297	volume_flow_l(1) = 0.0_wp
298
299	IF ( bc_radiation_l ) THEN
300	i = 0
301	ELSEIF ( bc_radiation_r ) THEN
302	i = nx+1
303	ENDIF
304
305	DO j = nys, nyn
306	!
307	!-- Sum up the volume flow through the south/north boundary
308	DO k = nzb+1, nzt
309	volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k) &
310	* MERGE( 1.0_wp, 0.0_wp, &
311	BTEST( wall_flags_0(k,j,i), 1 ) &
312	)
313	ENDDO
314	ENDDO
315
316	#if defined( __parallel )
317	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
318	CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 1, MPI_REAL, &
319	MPI_SUM, comm1dy, ierr )
320	#else
321	volume_flow = volume_flow_l
322	#endif
323	volume_flow_offset(1) = ( volume_flow_initial(1) - volume_flow(1) ) &
324	/ volume_flow_area(1)
325
326	DO j = nysg, nyng
327	DO k = nzb+1, nzt
328	u(k,j,i) = u(k,j,i) + volume_flow_offset(1) &
329	* MERGE( 1.0_wp, 0.0_wp, &
330	BTEST( wall_flags_0(k,j,i), 1 ) &
331	)
332	ENDDO
333	ENDDO
334
335	ENDIF
336
337	!
338	!-- South/north
339	IF ( conserve_volume_flow .AND. ( bc_radiation_n .OR. bc_radiation_s ) ) THEN
340
341	volume_flow(2) = 0.0_wp
342	volume_flow_l(2) = 0.0_wp
343
344	IF ( bc_radiation_s ) THEN
345	j = 0
346	ELSEIF ( bc_radiation_n ) THEN
347	j = ny+1
348	ENDIF
349
350	DO i = nxl, nxr
351	!
352	!-- Sum up the volume flow through the south/north boundary
353	DO k = nzb+1, nzt
354	volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k) &
355	* MERGE( 1.0_wp, 0.0_wp, &
356	BTEST( wall_flags_0(k,j,i), 2 ) &
357	)
358	ENDDO
359	ENDDO
360
361	#if defined( __parallel )
362	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
363	CALL MPI_ALLREDUCE( volume_flow_l(2), volume_flow(2), 1, MPI_REAL, &
364	MPI_SUM, comm1dx, ierr )
365	#else
366	volume_flow = volume_flow_l
367	#endif
368	volume_flow_offset(2) = ( volume_flow_initial(2) - volume_flow(2) ) &
369	/ volume_flow_area(2)
370
371	DO i = nxlg, nxrg
372	DO k = nzb+1, nzt
373	v(k,j,i) = v(k,j,i) + volume_flow_offset(2) &
374	* MERGE( 1.0_wp, 0.0_wp, &
375	BTEST( wall_flags_0(k,j,i), 2 ) &
376	)
377	ENDDO
378	ENDDO
379
380	ENDIF
381
382	!
383	!-- Remove mean vertical velocity in case that Neumann conditions are
384	!-- used both at bottom and top boundary, and if not a nested domain in a
385	!-- normal nesting run. In case of vertical nesting, this must be done.
386	!-- Therefore an auxiliary logical variable child_domain_nvn is used here, and
387	!-- nvn stands for non-vertical nesting.
388	!-- This cannot be done before the first initial time step because ngp_2dh_outer
389	!-- is not yet known then.
390	child_domain_nvn = child_domain
391	IF ( child_domain .AND. nesting_mode == 'vertical' ) THEN
392	child_domain_nvn = .FALSE.
393	ENDIF
394
395	IF ( ibc_p_b == 1 .AND. ibc_p_t == 1 .AND. .NOT. nesting_offline .AND. &
396	.NOT. child_domain_nvn .AND. intermediate_timestep_count /= 0 ) &
397	THEN
398	w_l = 0.0_wp; w_l_l = 0.0_wp
399	DO i = nxl, nxr
400	DO j = nys, nyn
401	DO k = nzb+1, nzt
402	w_l_l(k) = w_l_l(k) + w(k,j,i) &
403	* MERGE( 1.0_wp, 0.0_wp, &
404	BTEST( wall_flags_0(k,j,i), 3 ) &
405	)
406	ENDDO
407	ENDDO
408	ENDDO
409	#if defined( __parallel )
410	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
411	CALL MPI_ALLREDUCE( w_l_l(1), w_l(1), nzt, MPI_REAL, MPI_SUM, &
412	comm2d, ierr )
413	#else
414	w_l = w_l_l
415	#endif
416	DO k = 1, nzt
417	w_l(k) = w_l(k) / ngp_2dh_outer(k,0)
418	ENDDO
419	DO i = nxlg, nxrg
420	DO j = nysg, nyng
421	DO k = nzb+1, nzt
422	w(k,j,i) = w(k,j,i) - w_l(k) &
423	* MERGE( 1.0_wp, 0.0_wp, &
424	BTEST( wall_flags_0(k,j,i), 3 ) &
425	)
426	ENDDO
427	ENDDO
428	ENDDO
429	ENDIF
430
431	!
432	!-- Compute the divergence of the provisional velocity field.
433	CALL cpu_log( log_point_s(1), 'divergence', 'start' )
434
435	IF ( psolver(1:9) == 'multigrid' ) THEN
436	!$OMP PARALLEL DO SCHEDULE( STATIC ) PRIVATE (i,j,k)
437	DO i = nxl-1, nxr+1
438	DO j = nys-1, nyn+1
439	DO k = nzb, nzt+1
440	d(k,j,i) = 0.0_wp
441	ENDDO
442	ENDDO
443	ENDDO
444	ELSE
445	!$OMP PARALLEL DO SCHEDULE( STATIC ) PRIVATE (i,j,k)
446	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
447	!$ACC PRESENT(d)
448	DO i = nxl, nxr
449	DO j = nys, nyn
450	DO k = nzb+1, nzt
451	d(k,j,i) = 0.0_wp
452	ENDDO
453	ENDDO
454	ENDDO
455	ENDIF
456
457	localsum = 0.0_wp
458	threadsum = 0.0_wp
459
460	#if defined( __ibm )
461	!$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
462	!$OMP DO SCHEDULE( STATIC )
463	DO i = nxl, nxr
464	DO j = nys, nyn
465	DO k = nzb+1, nzt
466	d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx + &
467	( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy + &
468	( w(k,j,i) * rho_air_zw(k) - &
469	w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k) &
470	) * ddt_3d * d_weight_pres &
471	* MERGE( 1.0_wp, 0.0_wp, &
472	BTEST( wall_flags_0(k,j,i), 0 ) &
473	)
474	ENDDO
475	!
476	!-- Compute possible PE-sum of divergences for flow_statistics
477	DO k = nzb+1, nzt
478	threadsum = threadsum + ABS( d(k,j,i) ) &
479	* MERGE( 1.0_wp, 0.0_wp, &
480	BTEST( wall_flags_0(k,j,i), 0 ) &
481	)
482	ENDDO
483
484	ENDDO
485	ENDDO
486
487	IF ( intermediate_timestep_count == intermediate_timestep_count_max .OR. &
488	intermediate_timestep_count == 0 ) THEN
489	localsum = localsum + threadsum * dt_3d * weight_pres_l
490	ENDIF
491	!$OMP END PARALLEL
492	#else
493
494	!$OMP PARALLEL PRIVATE (i,j,k)
495	!$OMP DO SCHEDULE( STATIC )
496	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
497	!$ACC PRESENT(u, v, w, rho_air, rho_air_zw, ddzw, wall_flags_0) &
498	!$ACC PRESENT(d)
499	DO i = nxl, nxr
500	DO j = nys, nyn
501	DO k = 1, nzt
502	d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx + &
503	( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy + &
504	( w(k,j,i) * rho_air_zw(k) - &
505	w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k) &
506	) * ddt_3d * d_weight_pres &
507	* MERGE( 1.0_wp, 0.0_wp, &
508	BTEST( wall_flags_0(k,j,i), 0 ) &
509	)
510	ENDDO
511	ENDDO
512	ENDDO
513	!$OMP END PARALLEL
514
515	!
516	!-- Compute possible PE-sum of divergences for flow_statistics. Carry out
517	!-- computation only at last Runge-Kutta substep.
518	IF ( intermediate_timestep_count == intermediate_timestep_count_max .OR. &
519	intermediate_timestep_count == 0 ) THEN
520	!$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
521	!$OMP DO SCHEDULE( STATIC )
522	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
523	!$ACC REDUCTION(+:threadsum) COPY(threadsum) &
524	!$ACC PRESENT(d)
525	DO i = nxl, nxr
526	DO j = nys, nyn
527	DO k = nzb+1, nzt
528	threadsum = threadsum + ABS( d(k,j,i) )
529	ENDDO
530	ENDDO
531	ENDDO
532	localsum = localsum + threadsum * dt_3d * weight_pres_l
533	!$OMP END PARALLEL
534	ENDIF
535	#endif
536
537	!
538	!-- For completeness, set the divergence sum of all statistic regions to those
539	!-- of the total domain
540	IF ( intermediate_timestep_count == intermediate_timestep_count_max .OR. &
541	intermediate_timestep_count == 0 ) THEN
542	sums_divold_l(0:statistic_regions) = localsum
543	ENDIF
544
545	CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
546
547	!
548	!-- Compute the pressure perturbation solving the Poisson equation
549	IF ( psolver == 'poisfft' ) THEN
550
551	!
552	!-- Solve Poisson equation via FFT and solution of tridiagonal matrices
553	CALL poisfft( d )
554
555	!
556	!-- Store computed perturbation pressure and set boundary condition in
557	!-- z-direction
558	!$OMP PARALLEL DO PRIVATE (i,j,k)
559	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
560	!$ACC PRESENT(d, tend)
561	DO i = nxl, nxr
562	DO j = nys, nyn
563	DO k = nzb+1, nzt
564	tend(k,j,i) = d(k,j,i)
565	ENDDO
566	ENDDO
567	ENDDO
568
569	!
570	!-- Bottom boundary:
571	!-- This condition is only required for internal output. The pressure
572	!-- gradient (dp(nzb+1)-dp(nzb))/dz is not used anywhere else.
573	IF ( ibc_p_b == 1 ) THEN
574	!
575	!-- Neumann (dp/dz = 0). Using surfae data type, first for non-natural
576	!-- surfaces, then for natural and urban surfaces
577	!-- Upward facing
578	!$OMP PARALLEL DO PRIVATE( i, j, k )
579	!$ACC PARALLEL LOOP PRIVATE(i, j, k) &
580	!$ACC PRESENT(bc_h, tend)
581	DO m = 1, bc_h(0)%ns
582	i = bc_h(0)%i(m)
583	j = bc_h(0)%j(m)
584	k = bc_h(0)%k(m)
585	tend(k-1,j,i) = tend(k,j,i)
586	ENDDO
587	!
588	!-- Downward facing
589	!$OMP PARALLEL DO PRIVATE( i, j, k )
590	!$ACC PARALLEL LOOP PRIVATE(i, j, k) &
591	!$ACC PRESENT(bc_h, tend)
592	DO m = 1, bc_h(1)%ns
593	i = bc_h(1)%i(m)
594	j = bc_h(1)%j(m)
595	k = bc_h(1)%k(m)
596	tend(k+1,j,i) = tend(k,j,i)
597	ENDDO
598
599	ELSE
600	!
601	!-- Dirichlet. Using surface data type, first for non-natural
602	!-- surfaces, then for natural and urban surfaces
603	!-- Upward facing
604	!$OMP PARALLEL DO PRIVATE( i, j, k )
605	DO m = 1, bc_h(0)%ns
606	i = bc_h(0)%i(m)
607	j = bc_h(0)%j(m)
608	k = bc_h(0)%k(m)
609	tend(k-1,j,i) = 0.0_wp
610	ENDDO
611	!
612	!-- Downward facing
613	!$OMP PARALLEL DO PRIVATE( i, j, k )
614	DO m = 1, bc_h(1)%ns
615	i = bc_h(1)%i(m)
616	j = bc_h(1)%j(m)
617	k = bc_h(1)%k(m)
618	tend(k+1,j,i) = 0.0_wp
619	ENDDO
620
621	ENDIF
622
623	!
624	!-- Top boundary
625	IF ( ibc_p_t == 1 ) THEN
626	!
627	!-- Neumann
628	!$OMP PARALLEL DO PRIVATE (i,j,k)
629	DO i = nxlg, nxrg
630	DO j = nysg, nyng
631	tend(nzt+1,j,i) = tend(nzt,j,i)
632	ENDDO
633	ENDDO
634
635	ELSE
636	!
637	!-- Dirichlet
638	!$OMP PARALLEL DO PRIVATE (i,j,k)
639	!$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
640	!$ACC PRESENT(tend)
641	DO i = nxlg, nxrg
642	DO j = nysg, nyng
643	tend(nzt+1,j,i) = 0.0_wp
644	ENDDO
645	ENDDO
646
647	ENDIF
648
649	!
650	!-- Exchange boundaries for p
651	CALL exchange_horiz( tend, nbgp )
652
653	ELSEIF ( psolver == 'sor' ) THEN
654
655	!
656	!-- Solve Poisson equation for perturbation pressure using SOR-Red/Black
657	!-- scheme
658	CALL sor( d, ddzu_pres, ddzw, p_loc )
659	tend = p_loc
660
661	ELSEIF ( psolver(1:9) == 'multigrid' ) THEN
662
663	!
664	!-- Solve Poisson equation for perturbation pressure using Multigrid scheme,
665	!-- array tend is used to store the residuals.
666
667	!-- If the number of grid points of the gathered grid, which is collected
668	!-- on PE0, is larger than the number of grid points of an PE, than array
669	!-- tend will be enlarged.
670	IF ( gathered_size > subdomain_size ) THEN
671	DEALLOCATE( tend )
672	ALLOCATE( tend(nzb:nzt_mg(mg_switch_to_pe0_level)+1,nys_mg( &
673	mg_switch_to_pe0_level)-1:nyn_mg(mg_switch_to_pe0_level)+1,&
674	nxl_mg(mg_switch_to_pe0_level)-1:nxr_mg( &
675	mg_switch_to_pe0_level)+1) )
676	ENDIF
677
678	IF ( psolver == 'multigrid' ) THEN
679	CALL poismg( tend )
680	ELSE
681	CALL poismg_noopt( tend )
682	ENDIF
683
684	IF ( gathered_size > subdomain_size ) THEN
685	DEALLOCATE( tend )
686	ALLOCATE( tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
687	ENDIF
688
689	!
690	!-- Restore perturbation pressure on tend because this array is used
691	!-- further below to correct the velocity fields
692	DO i = nxl-1, nxr+1
693	DO j = nys-1, nyn+1
694	DO k = nzb, nzt+1
695	tend(k,j,i) = p_loc(k,j,i)
696	ENDDO
697	ENDDO
698	ENDDO
699
700	ENDIF
701
702	!
703	!-- Store perturbation pressure on array p, used for pressure data output.
704	!-- Ghost layers are added in the output routines (except sor-method: see below)
705	IF ( intermediate_timestep_count <= 1 ) THEN
706	!$OMP PARALLEL PRIVATE (i,j,k)
707	!$OMP DO
708	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
709	!$ACC PRESENT(p, tend)
710	DO i = nxl-1, nxr+1
711	DO j = nys-1, nyn+1
712	DO k = nzb, nzt+1
713	p(k,j,i) = tend(k,j,i) * &
714	weight_substep_l
715	ENDDO
716	ENDDO
717	ENDDO
718	!$OMP END PARALLEL
719
720	ELSEIF ( intermediate_timestep_count > 1 ) THEN
721	!$OMP PARALLEL PRIVATE (i,j,k)
722	!$OMP DO
723	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
724	!$ACC PRESENT(p, tend)
725	DO i = nxl-1, nxr+1
726	DO j = nys-1, nyn+1
727	DO k = nzb, nzt+1
728	p(k,j,i) = p(k,j,i) + tend(k,j,i) * &
729	weight_substep_l
730	ENDDO
731	ENDDO
732	ENDDO
733	!$OMP END PARALLEL
734
735	ENDIF
736
737	!
738	!-- SOR-method needs ghost layers for the next timestep
739	IF ( psolver == 'sor' ) CALL exchange_horiz( p, nbgp )
740
741	!
742	!-- Correction of the provisional velocities with the current perturbation
743	!-- pressure just computed
744	IF ( conserve_volume_flow .AND. ( bc_lr_cyc .OR. bc_ns_cyc ) ) THEN
745	volume_flow_l(1) = 0.0_wp
746	volume_flow_l(2) = 0.0_wp
747	ENDIF
748	!
749	!-- Add pressure gradients to the velocity components. Note, the loops are
750	!-- running over the entire model domain, even though, in case of non-cyclic
751	!-- boundaries u- and v-component are not prognostic at i=0 and j=0,
752	!-- respectiveley. However, in case of Dirichlet boundary conditions for the
753	!-- velocities, zero-gradient conditions for the pressure are set, so that
754	!-- no modification is imposed at the boundaries.
755	!$OMP PARALLEL PRIVATE (i,j,k)
756	!$OMP DO
757	!$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k) &
758	!$ACC PRESENT(u, v, w, tend, ddzu, wall_flags_0)
759	DO i = nxl, nxr
760	DO j = nys, nyn
761
762	DO k = nzb+1, nzt
763	w(k,j,i) = w(k,j,i) - dt_3d * &
764	( tend(k+1,j,i) - tend(k,j,i) ) * ddzu(k+1) &
765	* weight_pres_l &
766	* MERGE( 1.0_wp, 0.0_wp, &
767	BTEST( wall_flags_0(k,j,i), 3 ) &
768	)
769	ENDDO
770
771	DO k = nzb+1, nzt
772	u(k,j,i) = u(k,j,i) - dt_3d * &
773	( tend(k,j,i) - tend(k,j,i-1) ) * ddx &
774	* weight_pres_l &
775	* MERGE( 1.0_wp, 0.0_wp, &
776	BTEST( wall_flags_0(k,j,i), 1 ) &
777	)
778	ENDDO
779
780	DO k = nzb+1, nzt
781	v(k,j,i) = v(k,j,i) - dt_3d * &
782	( tend(k,j,i) - tend(k,j-1,i) ) * ddy &
783	* weight_pres_l &
784	* MERGE( 1.0_wp, 0.0_wp, &
785	BTEST( wall_flags_0(k,j,i), 2 ) &
786	)
787	ENDDO
788
789	ENDDO
790	ENDDO
791	!$OMP END PARALLEL
792
793	!
794	!-- The vertical velocity is not set to zero at nzt + 1 for nested domains
795	!-- Instead it is set to the values of nzt (see routine vnest_boundary_conds
796	!-- or pmci_interp_tril_t) BEFORE calling the pressure solver. To avoid jumps
797	!-- while plotting profiles w at the top has to be set to the values in the
798	!-- height nzt after above modifications. Hint: w level nzt+1 does not impact
799	!-- results.
800	IF ( child_domain .OR. coupling_mode == 'vnested_fine' ) THEN
801	w(nzt+1,:,:) = w(nzt,:,:)
802	ENDIF
803
804	!
805	!-- Sum up the volume flow through the right and north boundary
806	IF ( conserve_volume_flow .AND. bc_lr_cyc .AND. bc_ns_cyc .AND. &
807	nxr == nx ) THEN
808
809	!$OMP PARALLEL PRIVATE (j,k)
810	!$OMP DO
811	DO j = nys, nyn
812	!$OMP CRITICAL
813	DO k = nzb+1, nzt
814	volume_flow_l(1) = volume_flow_l(1) + u(k,j,nxr) * dzw(k) &
815	* MERGE( 1.0_wp, 0.0_wp, &
816	BTEST( wall_flags_0(k,j,nxr), 1 )&
817	)
818	ENDDO
819	!$OMP END CRITICAL
820	ENDDO
821	!$OMP END PARALLEL
822
823	ENDIF
824
825	IF ( conserve_volume_flow .AND. bc_ns_cyc .AND. bc_lr_cyc .AND. &
826	nyn == ny ) THEN
827
828	!$OMP PARALLEL PRIVATE (i,k)
829	!$OMP DO
830	DO i = nxl, nxr
831	!$OMP CRITICAL
832	DO k = nzb+1, nzt
833	volume_flow_l(2) = volume_flow_l(2) + v(k,nyn,i) * dzw(k) &
834	* MERGE( 1.0_wp, 0.0_wp, &
835	BTEST( wall_flags_0(k,nyn,i), 2 )&
836	)
837	ENDDO
838	!$OMP END CRITICAL
839	ENDDO
840	!$OMP END PARALLEL
841
842	ENDIF
843
844	!
845	!-- Conserve the volume flow
846	IF ( conserve_volume_flow .AND. ( bc_lr_cyc .AND. bc_ns_cyc ) ) THEN
847
848	#if defined( __parallel )
849	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
850	CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 2, MPI_REAL, &
851	MPI_SUM, comm2d, ierr )
852	#else
853	volume_flow = volume_flow_l
854	#endif
855
856	volume_flow_offset(1:2) = ( volume_flow_initial(1:2) - volume_flow(1:2) ) / &
857	volume_flow_area(1:2)
858
859	!$OMP PARALLEL PRIVATE (i,j,k)
860	!$OMP DO
861	DO i = nxl, nxr
862	DO j = nys, nyn
863	DO k = nzb+1, nzt
864	u(k,j,i) = u(k,j,i) + volume_flow_offset(1) &
865	* MERGE( 1.0_wp, 0.0_wp, &
866	BTEST( wall_flags_0(k,j,i), 1 ) &
867	)
868	ENDDO
869	DO k = nzb+1, nzt
870	v(k,j,i) = v(k,j,i) + volume_flow_offset(2) &
871	* MERGE( 1.0_wp, 0.0_wp, &
872	BTEST( wall_flags_0(k,j,i), 2 ) &
873	)
874	ENDDO
875	ENDDO
876	ENDDO
877
878	!$OMP END PARALLEL
879
880	ENDIF
881
882	!
883	!-- Exchange of boundaries for the velocities
884	CALL exchange_horiz( u, nbgp )
885	CALL exchange_horiz( v, nbgp )
886	CALL exchange_horiz( w, nbgp )
887
888	!
889	!-- Compute the divergence of the corrected velocity field,
890	!-- A possible PE-sum is computed in flow_statistics. Carry out computation
891	!-- only at last Runge-Kutta step.
892	IF ( intermediate_timestep_count == intermediate_timestep_count_max .OR. &
893	intermediate_timestep_count == 0 ) THEN
894	CALL cpu_log( log_point_s(1), 'divergence', 'start' )
895	sums_divnew_l = 0.0_wp
896
897	!
898	!-- d must be reset to zero because it can contain nonzero values below the
899	!-- topography
900	IF ( topography /= 'flat' ) d = 0.0_wp
901
902	localsum = 0.0_wp
903	threadsum = 0.0_wp
904
905	!$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
906	#if defined( __ibm )
907	!$OMP DO SCHEDULE( STATIC )
908	DO i = nxl, nxr
909	DO j = nys, nyn
910	DO k = nzb+1, nzt
911	d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx + &
912	( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy + &
913	( w(k,j,i) * rho_air_zw(k) - &
914	w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k) &
915	) * MERGE( 1.0_wp, 0.0_wp, &
916	BTEST( wall_flags_0(k,j,i), 0 ) &
917	)
918	ENDDO
919	DO k = nzb+1, nzt
920	threadsum = threadsum + ABS( d(k,j,i) )
921	ENDDO
922	ENDDO
923	ENDDO
924	#else
925	!$OMP DO SCHEDULE( STATIC )
926	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
927	!$ACC PRESENT(u, v, w, rho_air, rho_air_zw, ddzw, wall_flags_0) &
928	!$ACC PRESENT(d)
929	DO i = nxl, nxr
930	DO j = nys, nyn
931	DO k = nzb+1, nzt
932	d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx + &
933	( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy + &
934	( w(k,j,i) * rho_air_zw(k) - &
935	w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k) &
936	) * MERGE( 1.0_wp, 0.0_wp, &
937	BTEST( wall_flags_0(k,j,i), 0 ) &
938	)
939	ENDDO
940	ENDDO
941	ENDDO
942	!
943	!-- Compute possible PE-sum of divergences for flow_statistics
944	!$OMP DO SCHEDULE( STATIC )
945	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
946	!$ACC REDUCTION(+:threadsum) COPY(threadsum) &
947	!$ACC PRESENT(d)
948	DO i = nxl, nxr
949	DO j = nys, nyn
950	DO k = nzb+1, nzt
951	threadsum = threadsum + ABS( d(k,j,i) )
952	ENDDO
953	ENDDO
954	ENDDO
955	#endif
956
957	localsum = localsum + threadsum
958	!$OMP END PARALLEL
959
960	!
961	!-- For completeness, set the divergence sum of all statistic regions to those
962	!-- of the total domain
963	sums_divnew_l(0:statistic_regions) = localsum
964
965	CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
966
967	ENDIF
968
969	CALL cpu_log( log_point(8), 'pres', 'stop' )
970
971	END SUBROUTINE pres

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |