Home

Context Navigation

pres.f90 @ 92

Last change on this file since 92 was 90, checked in by raasch, 17 years ago

New:
---
Calculation and output of user-defined profiles. New &userpar parameters data_output_pr_user and max_pr_user.

check_parameters, flow_statistics, modules, parin, read_var_list, user_interface, write_var_list

Changed:

Division through dt_3d replaced by multiplication of the inverse. For performance optimisation, this is done in the loop calculating the divergence instead of using a seperate loop. (pres.f90) var_hom and var_sum renamed pr_palm.

data_output_profiles, flow_statistics, init_3d_model, modules, parin, pres, read_var_list, run_control, time_integration

Errors:

Bugfix: work_fft*_vec removed from some PRIVATE-declarations (poisfft).

Bugfix: field_chr renamed field_char (user_interface).

Bugfix: output of use_upstream_for_tke (header).

header, poisfft, user_interface

Property svn:keywords set to Id

File size: 17.3 KB

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

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

Context Navigation

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

Download in other formats: