Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |