Home

Context Navigation

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

Last change on this file since 719 was 719, checked in by gryschka, 13 years ago
Bugfix in volume flow control for double cyclic boundary conditions
Property svn:keywords set to `Id`
File size: 22.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |