Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |