Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |