Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |