Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |