Home

Context Navigation

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

Last change on this file since 4068 was 4015, checked in by raasch, 5 years ago
all reals changed to double precision in order to work with 32-bit working precision, otherwise calculated time intervals would mostly give zero; variable child_domain_nvn eliminated
Property svn:keywords set to `Id`
File size: 33.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |