Home

Context Navigation

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

Last change on this file since 1762 was 1762, checked in by hellstea, 8 years ago
Introduction of nested domain system
Property svn:keywords set to `Id`
File size: 24.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |