Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |