Home

Context Navigation

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

Last change on this file since 1682 was 1682, checked in by knoop, 9 years ago
Code annotations made doxygen readable
Property svn:keywords set to `Id`
File size: 24.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |