Home

Context Navigation

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

Last change on this file since 108 was 106, checked in by raasch, 17 years ago
preliminary update of bugfixes and extensions for non-cyclic BCs
Property svn:keywords set to `Id`
File size: 18.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |