Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |