Home

Context Navigation

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

Last change on this file since 1096 was 1093, checked in by raasch, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 23.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |