Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |