Home

Context Navigation

pres.f90 @ 622

Last change on this file since 622 was 622, checked in by raasch, 14 years ago

New:
---

Optional barriers included in order to speed up collective operations
MPI_ALLTOALL and MPI_ALLREDUCE. This feature is controlled with new initial
parameter collective_wait. Default is .FALSE, but .TRUE. on SGI-type
systems. (advec_particles, advec_s_bc, buoyancy, check_for_restart,
cpu_statistics, data_output_2d, data_output_ptseries, flow_statistics,
global_min_max, inflow_turbulence, init_3d_model, init_particles, init_pegrid,
init_slope, parin, pres, poismg, set_particle_attributes, timestep,
read_var_list, user_statistics, write_compressed, write_var_list)

Adjustments for Kyushu Univ. (lcrte, ibmku). Concerning hybrid
(MPI/openMP) runs, the number of openMP threads per MPI tasks can now
be given as an argument to mrun-option -O. (mbuild, mrun, subjob)

Changed:

Initialization of the module command changed for SGI-ICE/lcsgi (mbuild, subjob)

Errors:

Property svn:keywords set to Id

File size: 19.0 KB

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

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

Context Navigation

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

Download in other formats: