Home

Context Navigation

pres.f90 @ 707

Last change on this file since 707 was 707, checked in by raasch, 13 years ago

New:
---

In case of multigrid method, on coarse grid levels, gathered data are
identically processed on all PEs (before, on PE0 only), so that the subsequent
scattering of data is not neccessary any more. (modules, init_pegrid, poismg)

Changed:

Calculation of weighted average of p is now handled in the same way
regardless of the number of ghost layers (advection scheme). (pres)

multigrid and sor method are using p_loc for iterative
advancements of pressure. p_sub removed. (init_3d_model, modules, poismg, pres, sor)

bc_lr and bc_ns replaced by bc_lr_dirrad, bc_lr_raddir, bc_ns_dirrad, bc_ns_raddir
for speed optimization. (calc_spectra, check_parameters, exchange_horiz,
exchange_horiz_2d, header, init_3d_model, init_grid, init_pegrid, modules,
poismg, pres, sor, time_integration, timestep)

grid_level directly used as index for MPI data type arrays. (exchange_horiz,
poismg)

initial assignments of zero to array p for iterative solvers only (init_3d_model)

Errors:

localsum calculation modified for proper OpenMP reduction. (pres)

Bugfix: bottom (nzb) and top (nzt+1) boundary conditions set in routines
resid and restrict. They were missed before, which may have led to
unpredictable results. (poismg)

Property svn:keywords set to Id

File size: 22.2 KB

Line
1	SUBROUTINE pres
2
3	!------------------------------------------------------------------------------!
4	! Current revisions:
5	! -----------------
6	! Calculation of weighted average of p is now handled in the same way
7	! regardless of the number of ghost layers (advection scheme),
8	! multigrid and sor method are using p_loc for iterative advancements of
9	! pressure,
10	! localsum calculation modified for proper OpenMP reduction,
11	! bc_lr/ns replaced by bc_lr/ns_cyc
12	!
13	! Former revisions:
14	! -----------------
15	! $Id: pres.f90 707 2011-03-29 11:39:40Z raasch $
16	!
17	! 693 2011-03-08 09:..:..Z raasch
18	! bugfix: weighting coefficient added to ibm branch
19	!
20	! 680 2011-02-04 23:16:06Z gryschka
21	! bugfix: collective_wait
22	!
23	! 675 2011-01-19 10:56:55Z suehring
24	! Removed bugfix while copying tend.
25	!
26	! 673 2011-01-18 16:19:48Z suehring
27	! Weighting coefficients added for right computation of the pressure during
28	! Runge-Kutta substeps.
29	!
30	! 667 2010-12-23 12:06:00Z suehring/gryschka
31	! New allocation of tend when ws-scheme and multigrid is used. This is due to
32	! reasons of perforance of the data_exchange. The same is done with p after
33	! poismg is called.
34	! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng when no
35	! multigrid is used. Calls of exchange_horiz are modified.
36	! bugfix: After pressure correction no volume flow correction in case of
37	! non-cyclic boundary conditions
38	! (has to be done only before pressure correction)
39	! Call of SOR routine is referenced with ddzu_pres.
40	!
41	! 622 2010-12-10 08:08:13Z raasch
42	! optional barriers included in order to speed up collective operations
43	!
44	! 151 2008-03-07 13:42:18Z raasch
45	! Bugfix in volume flow control for non-cyclic boundary conditions
46	!
47	! 106 2007-08-16 14:30:26Z raasch
48	! Volume flow conservation added for the remaining three outflow boundaries
49	!
50	! 85 2007-05-11 09:35:14Z raasch
51	! Division through dt_3d replaced by multiplication of the inverse.
52	! For performance optimisation, this is done in the loop calculating the
53	! divergence instead of using a seperate loop.
54	!
55	! 75 2007-03-22 09:54:05Z raasch
56	! Volume flow control for non-cyclic boundary conditions added (currently only
57	! for the north boundary!!), 2nd+3rd argument removed from exchange horiz,
58	! mean vertical velocity is removed in case of Neumann boundary conditions
59	! both at the bottom and the top
60	!
61	! RCS Log replace by Id keyword, revision history cleaned up
62	!
63	! Revision 1.25 2006/04/26 13:26:12 raasch
64	! OpenMP optimization (+localsum, threadsum)
65	!
66	! Revision 1.1 1997/07/24 11:24:44 raasch
67	! Initial revision
68	!
69	!
70	! Description:
71	! ------------
72	! Compute the divergence of the provisional velocity field. Solve the Poisson
73	! equation for the perturbation pressure. Compute the final velocities using
74	! this perturbation pressure. Compute the remaining divergence.
75	!------------------------------------------------------------------------------!
76
77	USE arrays_3d
78	USE constants
79	USE control_parameters
80	USE cpulog
81	USE grid_variables
82	USE indices
83	USE interfaces
84	USE pegrid
85	USE poisfft_mod
86	USE poisfft_hybrid_mod
87	USE statistics
88
89	IMPLICIT NONE
90
91	INTEGER :: i, j, k, sr
92
93	REAL :: ddt_3d, localsum, threadsum, d_weight_pres
94
95	REAL, DIMENSION(1:2) :: volume_flow_l, volume_flow_offset
96	REAL, DIMENSION(1:nzt) :: w_l, w_l_l
97
98
99	CALL cpu_log( log_point(8), 'pres', 'start' )
100
101
102	ddt_3d = 1.0 / dt_3d
103	d_weight_pres = 1. / weight_pres(intermediate_timestep_count)
104
105	!
106	!-- Multigrid method expects array d to have one ghost layer.
107	!--
108	IF ( psolver == 'multigrid' ) THEN
109
110	DEALLOCATE( d )
111	ALLOCATE( d(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
112
113	!
114	!-- Since p is later used to hold the weighted average of the substeps, it
115	!-- cannot be used in the iterative solver. Therefore, its initial value is
116	!-- stored on p_loc, which is then iteratively advanced in every substep.
117	IF ( intermediate_timestep_count == 1 ) THEN
118	DO i = nxl-1, nxr+1
119	DO j = nys-1, nyn+1
120	DO k = nzb, nzt+1
121	p_loc(k,j,i) = p(k,j,i)
122	ENDDO
123	ENDDO
124	ENDDO
125	ENDIF
126
127	ELSEIF ( psolver == 'sor' .AND. intermediate_timestep_count == 1 ) THEN
128
129	!
130	!-- Since p is later used to hold the weighted average of the substeps, it
131	!-- cannot be used in the iterative solver. Therefore, its initial value is
132	!-- stored on p_loc, which is then iteratively advanced in every substep.
133	p_loc = p
134
135	ENDIF
136
137	!
138	!-- Conserve the volume flow at the outflow in case of non-cyclic lateral
139	!-- boundary conditions
140	!-- WARNING: so far, this conservation does not work at the left/south
141	!-- boundary if the topography at the inflow differs from that at the
142	!-- outflow! For this case, volume_flow_area needs adjustment!
143	!
144	!-- Left/right
145
146	IF ( conserve_volume_flow .AND. ( outflow_l .OR. outflow_r ) ) THEN
147
148	volume_flow(1) = 0.0
149	volume_flow_l(1) = 0.0
150
151	IF ( outflow_l ) THEN
152	i = 0
153	ELSEIF ( outflow_r ) THEN
154	i = nx+1
155	ENDIF
156
157	DO j = nys, nyn
158	!
159	!-- Sum up the volume flow through the south/north boundary
160	DO k = nzb_2d(j,i) + 1, nzt
161	volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k)
162	ENDDO
163	ENDDO
164
165	#if defined( __parallel )
166	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
167	CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 1, MPI_REAL, &
168	MPI_SUM, comm1dy, ierr )
169	#else
170	volume_flow = volume_flow_l
171	#endif
172	volume_flow_offset(1) = ( volume_flow_initial(1) - volume_flow(1) ) &
173	/ volume_flow_area(1)
174
175	DO j = nysg, nyng
176	DO k = nzb_2d(j,i) + 1, nzt
177	u(k,j,i) = u(k,j,i) + volume_flow_offset(1)
178	ENDDO
179	ENDDO
180
181	ENDIF
182
183	!
184	!-- South/north
185	IF ( conserve_volume_flow .AND. ( outflow_n .OR. outflow_s ) ) THEN
186
187	volume_flow(2) = 0.0
188	volume_flow_l(2) = 0.0
189
190	IF ( outflow_s ) THEN
191	j = 0
192	ELSEIF ( outflow_n ) THEN
193	j = ny+1
194	ENDIF
195
196	DO i = nxl, nxr
197	!
198	!-- Sum up the volume flow through the south/north boundary
199	DO k = nzb_2d(j,i) + 1, nzt
200	volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k)
201	ENDDO
202	ENDDO
203
204	#if defined( __parallel )
205	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
206	CALL MPI_ALLREDUCE( volume_flow_l(2), volume_flow(2), 1, MPI_REAL, &
207	MPI_SUM, comm1dx, ierr )
208	#else
209	volume_flow = volume_flow_l
210	#endif
211	volume_flow_offset(2) = ( volume_flow_initial(2) - volume_flow(2) ) &
212	/ volume_flow_area(2)
213
214	DO i = nxlg, nxrg
215	DO k = nzb_v_inner(j,i) + 1, nzt
216	v(k,j,i) = v(k,j,i) + volume_flow_offset(2)
217	ENDDO
218	ENDDO
219
220	ENDIF
221
222	!
223	!-- Remove mean vertical velocity
224	IF ( ibc_p_b == 1 .AND. ibc_p_t == 1 ) THEN
225	IF ( simulated_time > 0.0 ) THEN ! otherwise nzb_w_inner is not yet known
226	w_l = 0.0; w_l_l = 0.0
227	DO i = nxl, nxr
228	DO j = nys, nyn
229	DO k = nzb_w_inner(j,i)+1, nzt
230	w_l_l(k) = w_l_l(k) + w(k,j,i)
231	ENDDO
232	ENDDO
233	ENDDO
234	#if defined( __parallel )
235	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
236	CALL MPI_ALLREDUCE( w_l_l(1), w_l(1), nzt, MPI_REAL, MPI_SUM, comm2d, &
237	ierr )
238	#else
239	w_l = w_l_l
240	#endif
241	DO k = 1, nzt
242	w_l(k) = w_l(k) / ngp_2dh_outer(k,0)
243	ENDDO
244	DO i = nxlg, nxrg
245	DO j = nysg, nyng
246	DO k = nzb_w_inner(j,i)+1, nzt
247	w(k,j,i) = w(k,j,i) - w_l(k)
248	ENDDO
249	ENDDO
250	ENDDO
251	ENDIF
252	ENDIF
253
254	!
255	!-- Compute the divergence of the provisional velocity field.
256	CALL cpu_log( log_point_s(1), 'divergence', 'start' )
257
258	IF ( psolver == 'multigrid' ) THEN
259	!$OMP PARALLEL DO SCHEDULE( STATIC )
260	DO i = nxl-1, nxr+1
261	DO j = nys-1, nyn+1
262	DO k = nzb, nzt+1
263	d(k,j,i) = 0.0
264	ENDDO
265	ENDDO
266	ENDDO
267	ELSE
268	!$OMP PARALLEL DO SCHEDULE( STATIC )
269	DO i = nxl, nxra
270	DO j = nys, nyna
271	DO k = nzb+1, nzta
272	d(k,j,i) = 0.0
273	ENDDO
274	ENDDO
275	ENDDO
276	ENDIF
277
278	localsum = 0.0
279	threadsum = 0.0
280
281	#if defined( __ibm )
282	!$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
283	!$OMP DO SCHEDULE( STATIC )
284	DO i = nxl, nxr
285	DO j = nys, nyn
286	DO k = nzb_s_inner(j,i)+1, nzt
287	d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
288	( v(k,j+1,i) - v(k,j,i) ) * ddy + &
289	( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ) * ddt_3d &
290	* d_weight_pres
291	ENDDO
292	!
293	!-- Additional pressure boundary condition at the bottom boundary for
294	!-- inhomogeneous Prandtl layer heat fluxes and temperatures, respectively
295	!-- dp/dz = -(dtau13/dx + dtau23/dy) + g*pt'/pt0.
296	!-- This condition must not be applied at the start of a run, because then
297	!-- flow_statistics has not yet been called and thus sums = 0.
298	IF ( ibc_p_b == 2 .AND. sums(nzb+1,4) /= 0.0 ) THEN
299	k = nzb_s_inner(j,i)
300	d(k+1,j,i) = d(k+1,j,i) + ( &
301	( usws(j,i+1) - usws(j,i) ) * ddx &
302	+ ( vsws(j+1,i) - vsws(j,i) ) * ddy &
303	- g * ( pt(k+1,j,i) - sums(k+1,4) ) / &
304	sums(k+1,4) &
305	) * ddzw(k+1) * ddt_3d * d_weight_pres
306	ENDIF
307
308	!
309	!-- Compute possible PE-sum of divergences for flow_statistics
310	DO k = nzb_s_inner(j,i)+1, nzt
311	threadsum = threadsum + ABS( d(k,j,i) )
312	ENDDO
313
314	ENDDO
315	ENDDO
316
317	localsum = localsum + threadsum * dt_3d * &
318	weight_pres(intermediate_timestep_count)
319
320	!$OMP END PARALLEL
321	#else
322	IF ( ibc_p_b == 2 .AND. sums(nzb+1,4) /= 0.0 ) THEN
323	!$OMP PARALLEL PRIVATE (i,j,k)
324	!$OMP DO SCHEDULE( STATIC )
325	DO i = nxl, nxr
326	DO j = nys, nyn
327	DO k = nzb_s_inner(j,i)+1, nzt
328	d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
329	( v(k,j+1,i) - v(k,j,i) ) * ddy + &
330	( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ) * ddt_3d &
331	* d_weight_pres
332	ENDDO
333	ENDDO
334	!
335	!-- Additional pressure boundary condition at the bottom boundary for
336	!-- inhomogeneous Prandtl layer heat fluxes and temperatures, respectively
337	!-- dp/dz = -(dtau13/dx + dtau23/dy) + g*pt'/pt0.
338	!-- This condition must not be applied at the start of a run, because then
339	!-- flow_statistics has not yet been called and thus sums = 0.
340	DO j = nys, nyn
341	k = nzb_s_inner(j,i)
342	d(k+1,j,i) = d(k+1,j,i) + ( &
343	( usws(j,i+1) - usws(j,i) ) * ddx &
344	+ ( vsws(j+1,i) - vsws(j,i) ) * ddy &
345	- g * ( pt(k+1,j,i) - sums(k+1,4) ) / &
346	sums(k+1,4) &
347	) * ddzw(k+1) * ddt_3d &
348	* d_weight_pres
349	ENDDO
350	ENDDO
351	!$OMP END PARALLEL
352
353	ELSE
354
355	!$OMP PARALLEL PRIVATE (i,j,k)
356	!$OMP DO SCHEDULE( STATIC )
357	DO i = nxl, nxr
358	DO j = nys, nyn
359	DO k = nzb_s_inner(j,i)+1, nzt
360	d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
361	( v(k,j+1,i) - v(k,j,i) ) * ddy + &
362	( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ) * ddt_3d &
363	* d_weight_pres
364	ENDDO
365	ENDDO
366	ENDDO
367	!$OMP END PARALLEL
368
369	ENDIF
370
371	!
372	!-- Compute possible PE-sum of divergences for flow_statistics
373	!$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
374	!$OMP DO SCHEDULE( STATIC )
375	DO i = nxl, nxr
376	DO j = nys, nyn
377	DO k = nzb+1, nzt
378	threadsum = threadsum + ABS( d(k,j,i) )
379	ENDDO
380	ENDDO
381	ENDDO
382	localsum = localsum + threadsum * dt_3d * &
383	weight_pres(intermediate_timestep_count)
384	!$OMP END PARALLEL
385	#endif
386
387	!
388	!-- For completeness, set the divergence sum of all statistic regions to those
389	!-- of the total domain
390	sums_divold_l(0:statistic_regions) = localsum
391
392	CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
393
394	!
395	!-- Compute the pressure perturbation solving the Poisson equation
396	IF ( psolver(1:7) == 'poisfft' ) THEN
397
398	!
399	!-- Enlarge the size of tend, used as a working array for the transpositions
400	IF ( nxra > nxr .OR. nyna > nyn .OR. nza > nz ) THEN
401	DEALLOCATE( tend )
402	ALLOCATE( tend(1:nza,nys:nyna,nxl:nxra) )
403	ENDIF
404
405	!
406	!-- Solve Poisson equation via FFT and solution of tridiagonal matrices
407	IF ( psolver == 'poisfft' ) THEN
408	!
409	!-- Solver for 2d-decomposition
410	CALL poisfft( d, tend )
411	ELSEIF ( psolver == 'poisfft_hybrid' ) THEN
412	!
413	!-- Solver for 1d-decomposition (using MPI and OpenMP).
414	!-- The old hybrid-solver is still included here, as long as there
415	!-- are some optimization problems in poisfft
416	CALL poisfft_hybrid( d )
417	ENDIF
418
419	!
420	!-- Resize tend to its normal size
421	IF ( nxra > nxr .OR. nyna > nyn .OR. nza > nz ) THEN
422	DEALLOCATE( tend )
423	ALLOCATE( tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
424	ENDIF
425
426	!
427	!-- Store computed perturbation pressure and set boundary condition in
428	!-- z-direction
429	!$OMP PARALLEL DO
430	DO i = nxl, nxr
431	DO j = nys, nyn
432	DO k = nzb+1, nzt
433	tend(k,j,i) = d(k,j,i)
434	ENDDO
435	ENDDO
436	ENDDO
437
438	!
439	!-- Bottom boundary:
440	!-- This condition is only required for internal output. The pressure
441	!-- gradient (dp(nzb+1)-dp(nzb))/dz is not used anywhere else.
442	IF ( ibc_p_b == 1 ) THEN
443	!
444	!-- Neumann (dp/dz = 0)
445	!$OMP PARALLEL DO
446	DO i = nxlg, nxrg
447	DO j = nysg, nyng
448	tend(nzb_s_inner(j,i),j,i) = tend(nzb_s_inner(j,i)+1,j,i)
449	ENDDO
450	ENDDO
451
452	ELSEIF ( ibc_p_b == 2 ) THEN
453	!
454	!-- Neumann condition for inhomogeneous surfaces,
455	!-- here currently still in the form of a zero gradient. Actually
456	!-- dp/dz = -(dtau13/dx + dtau23/dy) + g*pt'/pt0 would have to be used for
457	!-- the computation (cf. above: computation of divergences).
458	!$OMP PARALLEL DO
459	DO i = nxlg, nxrg
460	DO j = nysg, nyng
461	tend(nzb_s_inner(j,i),j,i) = tend(nzb_s_inner(j,i)+1,j,i)
462	ENDDO
463	ENDDO
464
465	ELSE
466	!
467	!-- Dirichlet
468	!$OMP PARALLEL DO
469	DO i = nxlg, nxrg
470	DO j = nysg, nyng
471	tend(nzb_s_inner(j,i),j,i) = 0.0
472	ENDDO
473	ENDDO
474
475	ENDIF
476
477	!
478	!-- Top boundary
479	IF ( ibc_p_t == 1 ) THEN
480	!
481	!-- Neumann
482	!$OMP PARALLEL DO
483	DO i = nxlg, nxrg
484	DO j = nysg, nyng
485	tend(nzt+1,j,i) = tend(nzt,j,i)
486	ENDDO
487	ENDDO
488
489	ELSE
490	!
491	!-- Dirichlet
492	!$OMP PARALLEL DO
493	DO i = nxlg, nxrg
494	DO j = nysg, nyng
495	tend(nzt+1,j,i) = 0.0
496	ENDDO
497	ENDDO
498
499	ENDIF
500
501	!
502	!-- Exchange boundaries for p
503	CALL exchange_horiz( tend, nbgp )
504
505	ELSEIF ( psolver == 'sor' ) THEN
506
507	!
508	!-- Solve Poisson equation for perturbation pressure using SOR-Red/Black
509	!-- scheme
510	CALL sor( d, ddzu_pres, ddzw, p_loc )
511	tend = p_loc
512
513	ELSEIF ( psolver == 'multigrid' ) THEN
514
515	!
516	!-- Solve Poisson equation for perturbation pressure using Multigrid scheme,
517	!-- array tend is used to store the residuals, logical exchange_mg is used
518	!-- to discern data exchange in multigrid ( 1 ghostpoint ) and normal grid
519	!-- ( nbgp ghost points ).
520	CALL poismg( tend )
521
522	!
523	!-- Restore perturbation pressure on tend because this array is used
524	!-- further below to correct the velocity fields
525	DO i = nxl-1, nxr+1
526	DO j = nys-1, nyn+1
527	DO k = nzb, nzt+1
528	tend(k,j,i) = p_loc(k,j,i)
529	ENDDO
530	ENDDO
531	ENDDO
532
533	ENDIF
534
535	!
536	!-- Store perturbation pressure on array p, used for pressure data output.
537	!-- Ghost layers are added in the output routines (except sor-method: see below)
538	IF ( intermediate_timestep_count == 1 ) THEN
539	!$OMP PARALLEL PRIVATE (i,j,k)
540	!$OMP DO
541	DO i = nxl-1, nxr+1
542	DO j = nys-1, nyn+1
543	DO k = nzb, nzt+1
544	p(k,j,i) = tend(k,j,i) * &
545	weight_substep(intermediate_timestep_count)
546	ENDDO
547	ENDDO
548	ENDDO
549	!$OMP END PARALLEL
550
551	ELSE
552	!$OMP PARALLEL PRIVATE (i,j,k)
553	!$OMP DO
554	DO i = nxl-1, nxr+1
555	DO j = nys-1, nyn+1
556	DO k = nzb, nzt+1
557	p(k,j,i) = p(k,j,i) + tend(k,j,i) * &
558	weight_substep(intermediate_timestep_count)
559	ENDDO
560	ENDDO
561	ENDDO
562	!$OMP END PARALLEL
563
564	ENDIF
565
566	!
567	!-- SOR-method needs ghost layers for the next timestep
568	IF ( psolver == 'sor' ) CALL exchange_horiz( p, nbgp )
569
570	!
571	!-- Correction of the provisional velocities with the current perturbation
572	!-- pressure just computed
573	IF ( conserve_volume_flow .AND. ( bc_lr_cyc .OR. bc_ns_cyc ) ) THEN
574	volume_flow_l(1) = 0.0
575	volume_flow_l(2) = 0.0
576	ENDIF
577
578	!$OMP PARALLEL PRIVATE (i,j,k)
579	!$OMP DO
580	DO i = nxl, nxr
581	DO j = nys, nyn
582	DO k = nzb_w_inner(j,i)+1, nzt
583	w(k,j,i) = w(k,j,i) - dt_3d * &
584	( tend(k+1,j,i) - tend(k,j,i) ) * ddzu(k+1) * &
585	weight_pres(intermediate_timestep_count)
586	ENDDO
587	DO k = nzb_u_inner(j,i)+1, nzt
588	u(k,j,i) = u(k,j,i) - dt_3d * &
589	( tend(k,j,i) - tend(k,j,i-1) ) * ddx * &
590	weight_pres(intermediate_timestep_count)
591	ENDDO
592	DO k = nzb_v_inner(j,i)+1, nzt
593	v(k,j,i) = v(k,j,i) - dt_3d * &
594	( tend(k,j,i) - tend(k,j-1,i) ) * ddy * &
595	weight_pres(intermediate_timestep_count)
596	ENDDO
597	!
598	!-- Sum up the volume flow through the right and north boundary
599	IF ( conserve_volume_flow .AND. bc_lr_cyc .AND. bc_ns_cyc .AND. &
600	i == nx ) THEN
601	!$OMP CRITICAL
602	DO k = nzb_2d(j,i) + 1, nzt
603	volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k)
604	ENDDO
605	!$OMP END CRITICAL
606	ENDIF
607	IF ( conserve_volume_flow .AND. bc_ns_cyc .AND. bc_lr_cyc .AND. &
608	j == ny ) THEN
609	!$OMP CRITICAL
610	DO k = nzb_2d(j,i) + 1, nzt
611	volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k)
612	ENDDO
613	!$OMP END CRITICAL
614	ENDIF
615
616	ENDDO
617	ENDDO
618	!$OMP END PARALLEL
619
620	!
621	!-- Conserve the volume flow
622	IF ( conserve_volume_flow .AND. ( bc_lr_cyc .AND. bc_ns_cyc ) ) THEN
623
624	#if defined( __parallel )
625	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
626	CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 2, MPI_REAL, &
627	MPI_SUM, comm2d, ierr )
628	#else
629	volume_flow = volume_flow_l
630	#endif
631
632	volume_flow_offset = ( volume_flow_initial - volume_flow ) / &
633	volume_flow_area
634
635	!$OMP PARALLEL PRIVATE (i,j,k)
636	!$OMP DO
637	DO i = nxl, nxr
638	DO j = nys, nyn
639	DO k = nzb_u_inner(j,i) + 1, nzt
640	u(k,j,i) = u(k,j,i) + volume_flow_offset(1)
641	v(k,j,i) = v(k,j,i) + volume_flow_offset(2)
642	ENDDO
643	ENDDO
644	ENDDO
645
646	!$OMP END PARALLEL
647
648	ENDIF
649
650	!
651	!-- Exchange of boundaries for the velocities
652	CALL exchange_horiz( u, nbgp )
653	CALL exchange_horiz( v, nbgp )
654	CALL exchange_horiz( w, nbgp )
655
656	!
657	!-- Compute the divergence of the corrected velocity field,
658	!-- a possible PE-sum is computed in flow_statistics
659	CALL cpu_log( log_point_s(1), 'divergence', 'start' )
660	sums_divnew_l = 0.0
661
662	!
663	!-- d must be reset to zero because it can contain nonzero values below the
664	!-- topography
665	IF ( topography /= 'flat' ) d = 0.0
666
667	localsum = 0.0
668	threadsum = 0.0
669
670	!$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
671	!$OMP DO SCHEDULE( STATIC )
672	#if defined( __ibm )
673	DO i = nxl, nxr
674	DO j = nys, nyn
675	DO k = nzb_s_inner(j,i)+1, nzt
676	d(k,j,i) = ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
677	( v(k,j+1,i) - v(k,j,i) ) * ddy + &
678	( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
679	ENDDO
680	DO k = nzb+1, nzt
681	threadsum = threadsum + ABS( d(k,j,i) )
682	ENDDO
683	ENDDO
684	ENDDO
685	#else
686	DO i = nxl, nxr
687	DO j = nys, nyn
688	DO k = nzb_s_inner(j,i)+1, nzt
689	d(k,j,i) = ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
690	( v(k,j+1,i) - v(k,j,i) ) * ddy + &
691	( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
692	threadsum = threadsum + ABS( d(k,j,i) )
693	ENDDO
694	ENDDO
695	ENDDO
696	#endif
697
698	localsum = localsum + threadsum
699	!$OMP END PARALLEL
700
701	!
702	!-- For completeness, set the divergence sum of all statistic regions to those
703	!-- of the total domain
704	sums_divnew_l(0:statistic_regions) = localsum
705
706	CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
707
708	CALL cpu_log( log_point(8), 'pres', 'stop' )
709
710
711
712	END SUBROUTINE pres

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

Context Navigation

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

Download in other formats: