Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |