Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |