Home

Context Navigation

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

Last change on this file since 680 was 680, checked in by gryschka, 13 years ago
message string
Property svn:keywords set to `Id`
File size: 23.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |