Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |