Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |