Home

Context Navigation

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

Last change on this file since 2250 was 2233, checked in by suehring, 7 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 30.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |