Home

Context Navigation

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

Last change on this file since 4566 was 4457, checked in by raasch, 5 years ago
ghost point exchange modularized, bugfix for wrong 2d-exchange
Property svn:keywords set to `Id`
File size: 30.3 KB

Line
1	!> @file pres.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of the PALM model system.
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-2020 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! ------------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: pres.f90 4457 2020-03-11 14:20:43Z suehring $
27	! use statement for exchange horiz added
28	!
29	! 4360 2020-01-07 11:25:50Z suehring
30	! Introduction of wall_flags_total_0, which currently sets bits based on static
31	! topography information used in wall_flags_static_0
32	!
33	! 4329 2019-12-10 15:46:36Z motisi
34	! Renamed wall_flags_0 to wall_flags_static_0
35	!
36	! 4182 2019-08-22 15:20:23Z scharf
37	! Corrected "Former revisions" section
38	!
39	! 4015 2019-06-05 13:25:35Z raasch
40	! variable child_domain_nvn eliminated
41	!
42	! 3849 2019-04-01 16:35:16Z knoop
43	! OpenACC port for SPEC
44	!
45	! Revision 1.1 1997/07/24 11:24:44 raasch
46	! Initial revision
47	!
48	!
49	! Description:
50	! ------------
51	!> Compute the divergence of the provisional velocity field. Solve the Poisson
52	!> equation for the perturbation pressure. Compute the final velocities using
53	!> this perturbation pressure. Compute the remaining divergence.
54	!------------------------------------------------------------------------------!
55	SUBROUTINE pres
56
57
58	USE arrays_3d, &
59	ONLY: d, ddzu, ddzu_pres, ddzw, dzw, p, p_loc, rho_air, rho_air_zw, &
60	tend, u, v, w
61
62	USE control_parameters, &
63	ONLY: bc_lr_cyc, bc_ns_cyc, bc_radiation_l, bc_radiation_n, &
64	bc_radiation_r, bc_radiation_s, child_domain, &
65	conserve_volume_flow, coupling_mode, &
66	dt_3d, gathered_size, ibc_p_b, ibc_p_t, &
67	intermediate_timestep_count, intermediate_timestep_count_max, &
68	mg_switch_to_pe0_level, nesting_offline, &
69	psolver, subdomain_size, &
70	topography, volume_flow, volume_flow_area, volume_flow_initial
71
72	USE cpulog, &
73	ONLY: cpu_log, log_point, log_point_s
74
75	USE exchange_horiz_mod, &
76	ONLY: exchange_horiz
77
78	USE grid_variables, &
79	ONLY: ddx, ddy
80
81	USE indices, &
82	ONLY: nbgp, ngp_2dh_outer, nx, nxl, nxlg, nxl_mg, nxr, nxrg, nxr_mg, &
83	ny, nys, nysg, nys_mg, nyn, nyng, nyn_mg, nzb, nzt, nzt_mg, &
84	wall_flags_total_0
85
86	USE kinds
87
88	USE pegrid
89
90	USE pmc_interface, &
91	ONLY: nesting_mode
92
93	USE poisfft_mod, &
94	ONLY: poisfft
95
96	USE poismg_mod
97
98	USE poismg_noopt_mod
99
100	USE statistics, &
101	ONLY: statistic_regions, sums_divnew_l, sums_divold_l, weight_pres, &
102	weight_substep
103
104	USE surface_mod, &
105	ONLY : bc_h
106
107	IMPLICIT NONE
108
109	INTEGER(iwp) :: i !<
110	INTEGER(iwp) :: j !<
111	INTEGER(iwp) :: k !<
112	INTEGER(iwp) :: m !<
113
114	REAL(wp) :: ddt_3d !<
115	REAL(wp) :: d_weight_pres !<
116	REAL(wp) :: localsum !<
117	REAL(wp) :: threadsum !<
118	REAL(wp) :: weight_pres_l !<
119	REAL(wp) :: weight_substep_l !<
120
121	REAL(wp), DIMENSION(1:3) :: volume_flow_l !<
122	REAL(wp), DIMENSION(1:3) :: volume_flow_offset !<
123	REAL(wp), DIMENSION(1:nzt) :: w_l !<
124	REAL(wp), DIMENSION(1:nzt) :: w_l_l !<
125
126
127	CALL cpu_log( log_point(8), 'pres', 'start' )
128
129	!
130	!-- Calculate quantities to be used locally
131	ddt_3d = 1.0_wp / dt_3d
132	IF ( intermediate_timestep_count == 0 ) THEN
133	!
134	!-- If pres is called before initial time step
135	weight_pres_l = 1.0_wp
136	d_weight_pres = 1.0_wp
137	weight_substep_l = 1.0_wp
138	ELSE
139	weight_pres_l = weight_pres(intermediate_timestep_count)
140	d_weight_pres = 1.0_wp / weight_pres(intermediate_timestep_count)
141	weight_substep_l = weight_substep(intermediate_timestep_count)
142	ENDIF
143
144	!
145	!-- Multigrid method expects array d to have one ghost layer.
146	!--
147	IF ( psolver(1:9) == 'multigrid' ) THEN
148
149	DEALLOCATE( d )
150	ALLOCATE( d(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
151
152	!
153	!-- Since p is later used to hold the weighted average of the substeps, it
154	!-- cannot be used in the iterative solver. Therefore, its initial value is
155	!-- stored on p_loc, which is then iteratively advanced in every substep.
156	IF ( intermediate_timestep_count <= 1 ) THEN
157	DO i = nxl-1, nxr+1
158	DO j = nys-1, nyn+1
159	DO k = nzb, nzt+1
160	p_loc(k,j,i) = p(k,j,i)
161	ENDDO
162	ENDDO
163	ENDDO
164	ENDIF
165
166	ELSEIF ( psolver == 'sor' .AND. intermediate_timestep_count <= 1 ) THEN
167
168	!
169	!-- Since p is later used to hold the weighted average of the substeps, it
170	!-- cannot be used in the iterative solver. Therefore, its initial value is
171	!-- stored on p_loc, which is then iteratively advanced in every substep.
172	p_loc = p
173
174	ENDIF
175
176	!
177	!-- Conserve the volume flow at the outflow in case of non-cyclic lateral
178	!-- boundary conditions
179	!-- WARNING: so far, this conservation does not work at the left/south
180	!-- boundary if the topography at the inflow differs from that at the
181	!-- outflow! For this case, volume_flow_area needs adjustment!
182	!
183	!-- Left/right
184	IF ( conserve_volume_flow .AND. ( bc_radiation_l .OR. &
185	bc_radiation_r ) ) THEN
186
187	volume_flow(1) = 0.0_wp
188	volume_flow_l(1) = 0.0_wp
189
190	IF ( bc_radiation_l ) THEN
191	i = 0
192	ELSEIF ( bc_radiation_r ) THEN
193	i = nx+1
194	ENDIF
195
196	DO j = nys, nyn
197	!
198	!-- Sum up the volume flow through the south/north boundary
199	DO k = nzb+1, nzt
200	volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k) &
201	* MERGE( 1.0_wp, 0.0_wp, &
202	BTEST( wall_flags_total_0(k,j,i), 1 ) &
203	)
204	ENDDO
205	ENDDO
206
207	#if defined( __parallel )
208	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
209	CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 1, MPI_REAL, &
210	MPI_SUM, comm1dy, ierr )
211	#else
212	volume_flow = volume_flow_l
213	#endif
214	volume_flow_offset(1) = ( volume_flow_initial(1) - volume_flow(1) ) &
215	/ volume_flow_area(1)
216
217	DO j = nysg, nyng
218	DO k = nzb+1, nzt
219	u(k,j,i) = u(k,j,i) + volume_flow_offset(1) &
220	* MERGE( 1.0_wp, 0.0_wp, &
221	BTEST( wall_flags_total_0(k,j,i), 1 ) &
222	)
223	ENDDO
224	ENDDO
225
226	ENDIF
227
228	!
229	!-- South/north
230	IF ( conserve_volume_flow .AND. ( bc_radiation_n .OR. bc_radiation_s ) ) THEN
231
232	volume_flow(2) = 0.0_wp
233	volume_flow_l(2) = 0.0_wp
234
235	IF ( bc_radiation_s ) THEN
236	j = 0
237	ELSEIF ( bc_radiation_n ) THEN
238	j = ny+1
239	ENDIF
240
241	DO i = nxl, nxr
242	!
243	!-- Sum up the volume flow through the south/north boundary
244	DO k = nzb+1, nzt
245	volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k) &
246	* MERGE( 1.0_wp, 0.0_wp, &
247	BTEST( wall_flags_total_0(k,j,i), 2 ) &
248	)
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+1, nzt
264	v(k,j,i) = v(k,j,i) + volume_flow_offset(2) &
265	* MERGE( 1.0_wp, 0.0_wp, &
266	BTEST( wall_flags_total_0(k,j,i), 2 ) &
267	)
268	ENDDO
269	ENDDO
270
271	ENDIF
272
273	!
274	!-- Remove mean vertical velocity in case that Neumann conditions are used both at bottom and top
275	!-- boundary. With Neumann conditions at both vertical boundaries, the solver cannot remove
276	!-- mean vertical velocities. They should be removed, because incompressibility requires that
277	!-- the vertical gradient of vertical velocity is zero. Since w=0 at the solid surface, it must be
278	!-- zero everywhere.
279	!-- This must not be done in case of a 3d-nesting child domain, because a mean vertical velocity
280	!-- can physically exist in such a domain.
281	!-- Also in case of offline nesting, mean vertical velocities may exist (and must not be removed),
282	!-- caused by horizontal divergence/convergence of the large scale flow that is prescribed at the
283	!-- side boundaries.
284	!-- The removal cannot be done before the first initial time step because ngp_2dh_outer
285	!-- is not yet known then.
286	IF ( ibc_p_b == 1 .AND. ibc_p_t == 1 .AND. .NOT. nesting_offline &
287	.AND. .NOT. ( child_domain .AND. nesting_mode /= 'vertical' ) &
288	.AND. intermediate_timestep_count /= 0 ) &
289	THEN
290	w_l = 0.0_wp; w_l_l = 0.0_wp
291	DO i = nxl, nxr
292	DO j = nys, nyn
293	DO k = nzb+1, nzt
294	w_l_l(k) = w_l_l(k) + w(k,j,i) &
295	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
296	ENDDO
297	ENDDO
298	ENDDO
299	#if defined( __parallel )
300	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
301	CALL MPI_ALLREDUCE( w_l_l(1), w_l(1), nzt, MPI_REAL, MPI_SUM, comm2d, ierr )
302	#else
303	w_l = w_l_l
304	#endif
305	DO k = 1, nzt
306	w_l(k) = w_l(k) / ngp_2dh_outer(k,0)
307	ENDDO
308	DO i = nxlg, nxrg
309	DO j = nysg, nyng
310	DO k = nzb+1, nzt
311	w(k,j,i) = w(k,j,i) - w_l(k) &
312	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
313	ENDDO
314	ENDDO
315	ENDDO
316	ENDIF
317
318	!
319	!-- Compute the divergence of the provisional velocity field.
320	CALL cpu_log( log_point_s(1), 'divergence', 'start' )
321
322	IF ( psolver(1:9) == 'multigrid' ) THEN
323	!$OMP PARALLEL DO SCHEDULE( STATIC ) PRIVATE (i,j,k)
324	DO i = nxl-1, nxr+1
325	DO j = nys-1, nyn+1
326	DO k = nzb, nzt+1
327	d(k,j,i) = 0.0_wp
328	ENDDO
329	ENDDO
330	ENDDO
331	ELSE
332	!$OMP PARALLEL DO SCHEDULE( STATIC ) PRIVATE (i,j,k)
333	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
334	!$ACC PRESENT(d)
335	DO i = nxl, nxr
336	DO j = nys, nyn
337	DO k = nzb+1, nzt
338	d(k,j,i) = 0.0_wp
339	ENDDO
340	ENDDO
341	ENDDO
342	ENDIF
343
344	localsum = 0.0_wp
345	threadsum = 0.0_wp
346
347	#if defined( __ibm )
348	!$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
349	!$OMP DO SCHEDULE( STATIC )
350	DO i = nxl, nxr
351	DO j = nys, nyn
352	DO k = nzb+1, nzt
353	d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx + &
354	( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy + &
355	( w(k,j,i) * rho_air_zw(k) - &
356	w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k) &
357	) * ddt_3d * d_weight_pres &
358	* MERGE( 1.0_wp, 0.0_wp, &
359	BTEST( wall_flags_total_0(k,j,i), 0 ) &
360	)
361	ENDDO
362	!
363	!-- Compute possible PE-sum of divergences for flow_statistics
364	DO k = nzb+1, nzt
365	threadsum = threadsum + ABS( d(k,j,i) ) &
366	* MERGE( 1.0_wp, 0.0_wp, &
367	BTEST( wall_flags_total_0(k,j,i), 0 ) &
368	)
369	ENDDO
370
371	ENDDO
372	ENDDO
373
374	IF ( intermediate_timestep_count == intermediate_timestep_count_max .OR. &
375	intermediate_timestep_count == 0 ) THEN
376	localsum = localsum + threadsum * dt_3d * weight_pres_l
377	ENDIF
378	!$OMP END PARALLEL
379	#else
380
381	!$OMP PARALLEL PRIVATE (i,j,k)
382	!$OMP DO SCHEDULE( STATIC )
383	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
384	!$ACC PRESENT(u, v, w, rho_air, rho_air_zw, ddzw, wall_flags_total_0) &
385	!$ACC PRESENT(d)
386	DO i = nxl, nxr
387	DO j = nys, nyn
388	DO k = 1, nzt
389	d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx + &
390	( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy + &
391	( w(k,j,i) * rho_air_zw(k) - &
392	w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k) &
393	) * ddt_3d * d_weight_pres &
394	* MERGE( 1.0_wp, 0.0_wp, &
395	BTEST( wall_flags_total_0(k,j,i), 0 ) &
396	)
397	ENDDO
398	ENDDO
399	ENDDO
400	!$OMP END PARALLEL
401
402	!
403	!-- Compute possible PE-sum of divergences for flow_statistics. Carry out
404	!-- computation only at last Runge-Kutta substep.
405	IF ( intermediate_timestep_count == intermediate_timestep_count_max .OR. &
406	intermediate_timestep_count == 0 ) THEN
407	!$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
408	!$OMP DO SCHEDULE( STATIC )
409	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
410	!$ACC REDUCTION(+:threadsum) COPY(threadsum) &
411	!$ACC PRESENT(d)
412	DO i = nxl, nxr
413	DO j = nys, nyn
414	DO k = nzb+1, nzt
415	threadsum = threadsum + ABS( d(k,j,i) )
416	ENDDO
417	ENDDO
418	ENDDO
419	localsum = localsum + threadsum * dt_3d * weight_pres_l
420	!$OMP END PARALLEL
421	ENDIF
422	#endif
423
424	!
425	!-- For completeness, set the divergence sum of all statistic regions to those
426	!-- of the total domain
427	IF ( intermediate_timestep_count == intermediate_timestep_count_max .OR. &
428	intermediate_timestep_count == 0 ) THEN
429	sums_divold_l(0:statistic_regions) = localsum
430	ENDIF
431
432	CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
433
434	!
435	!-- Compute the pressure perturbation solving the Poisson equation
436	IF ( psolver == 'poisfft' ) THEN
437
438	!
439	!-- Solve Poisson equation via FFT and solution of tridiagonal matrices
440	CALL poisfft( d )
441
442	!
443	!-- Store computed perturbation pressure and set boundary condition in
444	!-- z-direction
445	!$OMP PARALLEL DO PRIVATE (i,j,k)
446	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
447	!$ACC PRESENT(d, tend)
448	DO i = nxl, nxr
449	DO j = nys, nyn
450	DO k = nzb+1, nzt
451	tend(k,j,i) = d(k,j,i)
452	ENDDO
453	ENDDO
454	ENDDO
455
456	!
457	!-- Bottom boundary:
458	!-- This condition is only required for internal output. The pressure
459	!-- gradient (dp(nzb+1)-dp(nzb))/dz is not used anywhere else.
460	IF ( ibc_p_b == 1 ) THEN
461	!
462	!-- Neumann (dp/dz = 0). Using surfae data type, first for non-natural
463	!-- surfaces, then for natural and urban surfaces
464	!-- Upward facing
465	!$OMP PARALLEL DO PRIVATE( i, j, k )
466	!$ACC PARALLEL LOOP PRIVATE(i, j, k) &
467	!$ACC PRESENT(bc_h, tend)
468	DO m = 1, bc_h(0)%ns
469	i = bc_h(0)%i(m)
470	j = bc_h(0)%j(m)
471	k = bc_h(0)%k(m)
472	tend(k-1,j,i) = tend(k,j,i)
473	ENDDO
474	!
475	!-- Downward facing
476	!$OMP PARALLEL DO PRIVATE( i, j, k )
477	!$ACC PARALLEL LOOP PRIVATE(i, j, k) &
478	!$ACC PRESENT(bc_h, tend)
479	DO m = 1, bc_h(1)%ns
480	i = bc_h(1)%i(m)
481	j = bc_h(1)%j(m)
482	k = bc_h(1)%k(m)
483	tend(k+1,j,i) = tend(k,j,i)
484	ENDDO
485
486	ELSE
487	!
488	!-- Dirichlet. Using surface data type, first for non-natural
489	!-- surfaces, then for natural and urban surfaces
490	!-- Upward facing
491	!$OMP PARALLEL DO PRIVATE( i, j, k )
492	DO m = 1, bc_h(0)%ns
493	i = bc_h(0)%i(m)
494	j = bc_h(0)%j(m)
495	k = bc_h(0)%k(m)
496	tend(k-1,j,i) = 0.0_wp
497	ENDDO
498	!
499	!-- Downward facing
500	!$OMP PARALLEL DO PRIVATE( i, j, k )
501	DO m = 1, bc_h(1)%ns
502	i = bc_h(1)%i(m)
503	j = bc_h(1)%j(m)
504	k = bc_h(1)%k(m)
505	tend(k+1,j,i) = 0.0_wp
506	ENDDO
507
508	ENDIF
509
510	!
511	!-- Top boundary
512	IF ( ibc_p_t == 1 ) THEN
513	!
514	!-- Neumann
515	!$OMP PARALLEL DO PRIVATE (i,j,k)
516	DO i = nxlg, nxrg
517	DO j = nysg, nyng
518	tend(nzt+1,j,i) = tend(nzt,j,i)
519	ENDDO
520	ENDDO
521
522	ELSE
523	!
524	!-- Dirichlet
525	!$OMP PARALLEL DO PRIVATE (i,j,k)
526	!$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
527	!$ACC PRESENT(tend)
528	DO i = nxlg, nxrg
529	DO j = nysg, nyng
530	tend(nzt+1,j,i) = 0.0_wp
531	ENDDO
532	ENDDO
533
534	ENDIF
535
536	!
537	!-- Exchange boundaries for p
538	CALL exchange_horiz( tend, nbgp )
539
540	ELSEIF ( psolver == 'sor' ) THEN
541
542	!
543	!-- Solve Poisson equation for perturbation pressure using SOR-Red/Black
544	!-- scheme
545	CALL sor( d, ddzu_pres, ddzw, p_loc )
546	tend = p_loc
547
548	ELSEIF ( psolver(1:9) == 'multigrid' ) THEN
549
550	!
551	!-- Solve Poisson equation for perturbation pressure using Multigrid scheme,
552	!-- array tend is used to store the residuals.
553
554	!-- If the number of grid points of the gathered grid, which is collected
555	!-- on PE0, is larger than the number of grid points of an PE, than array
556	!-- tend will be enlarged.
557	IF ( gathered_size > subdomain_size ) THEN
558	DEALLOCATE( tend )
559	ALLOCATE( tend(nzb:nzt_mg(mg_switch_to_pe0_level)+1,nys_mg( &
560	mg_switch_to_pe0_level)-1:nyn_mg(mg_switch_to_pe0_level)+1,&
561	nxl_mg(mg_switch_to_pe0_level)-1:nxr_mg( &
562	mg_switch_to_pe0_level)+1) )
563	ENDIF
564
565	IF ( psolver == 'multigrid' ) THEN
566	CALL poismg( tend )
567	ELSE
568	CALL poismg_noopt( tend )
569	ENDIF
570
571	IF ( gathered_size > subdomain_size ) THEN
572	DEALLOCATE( tend )
573	ALLOCATE( tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
574	ENDIF
575
576	!
577	!-- Restore perturbation pressure on tend because this array is used
578	!-- further below to correct the velocity fields
579	DO i = nxl-1, nxr+1
580	DO j = nys-1, nyn+1
581	DO k = nzb, nzt+1
582	tend(k,j,i) = p_loc(k,j,i)
583	ENDDO
584	ENDDO
585	ENDDO
586
587	ENDIF
588
589	!
590	!-- Store perturbation pressure on array p, used for pressure data output.
591	!-- Ghost layers are added in the output routines (except sor-method: see below)
592	IF ( intermediate_timestep_count <= 1 ) THEN
593	!$OMP PARALLEL PRIVATE (i,j,k)
594	!$OMP DO
595	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
596	!$ACC PRESENT(p, tend)
597	DO i = nxl-1, nxr+1
598	DO j = nys-1, nyn+1
599	DO k = nzb, nzt+1
600	p(k,j,i) = tend(k,j,i) * &
601	weight_substep_l
602	ENDDO
603	ENDDO
604	ENDDO
605	!$OMP END PARALLEL
606
607	ELSEIF ( intermediate_timestep_count > 1 ) THEN
608	!$OMP PARALLEL PRIVATE (i,j,k)
609	!$OMP DO
610	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
611	!$ACC PRESENT(p, tend)
612	DO i = nxl-1, nxr+1
613	DO j = nys-1, nyn+1
614	DO k = nzb, nzt+1
615	p(k,j,i) = p(k,j,i) + tend(k,j,i) * &
616	weight_substep_l
617	ENDDO
618	ENDDO
619	ENDDO
620	!$OMP END PARALLEL
621
622	ENDIF
623
624	!
625	!-- SOR-method needs ghost layers for the next timestep
626	IF ( psolver == 'sor' ) CALL exchange_horiz( p, nbgp )
627
628	!
629	!-- Correction of the provisional velocities with the current perturbation
630	!-- pressure just computed
631	IF ( conserve_volume_flow .AND. ( bc_lr_cyc .OR. bc_ns_cyc ) ) THEN
632	volume_flow_l(1) = 0.0_wp
633	volume_flow_l(2) = 0.0_wp
634	ENDIF
635	!
636	!-- Add pressure gradients to the velocity components. Note, the loops are
637	!-- running over the entire model domain, even though, in case of non-cyclic
638	!-- boundaries u- and v-component are not prognostic at i=0 and j=0,
639	!-- respectiveley. However, in case of Dirichlet boundary conditions for the
640	!-- velocities, zero-gradient conditions for the pressure are set, so that
641	!-- no modification is imposed at the boundaries.
642	!$OMP PARALLEL PRIVATE (i,j,k)
643	!$OMP DO
644	!$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k) &
645	!$ACC PRESENT(u, v, w, tend, ddzu, wall_flags_total_0)
646	DO i = nxl, nxr
647	DO j = nys, nyn
648
649	DO k = nzb+1, nzt
650	w(k,j,i) = w(k,j,i) - dt_3d * &
651	( tend(k+1,j,i) - tend(k,j,i) ) * ddzu(k+1) &
652	* weight_pres_l &
653	* MERGE( 1.0_wp, 0.0_wp, &
654	BTEST( wall_flags_total_0(k,j,i), 3 ) &
655	)
656	ENDDO
657
658	DO k = nzb+1, nzt
659	u(k,j,i) = u(k,j,i) - dt_3d * &
660	( tend(k,j,i) - tend(k,j,i-1) ) * ddx &
661	* weight_pres_l &
662	* MERGE( 1.0_wp, 0.0_wp, &
663	BTEST( wall_flags_total_0(k,j,i), 1 ) &
664	)
665	ENDDO
666
667	DO k = nzb+1, nzt
668	v(k,j,i) = v(k,j,i) - dt_3d * &
669	( tend(k,j,i) - tend(k,j-1,i) ) * ddy &
670	* weight_pres_l &
671	* MERGE( 1.0_wp, 0.0_wp, &
672	BTEST( wall_flags_total_0(k,j,i), 2 ) &
673	)
674	ENDDO
675
676	ENDDO
677	ENDDO
678	!$OMP END PARALLEL
679
680	!
681	!-- The vertical velocity is not set to zero at nzt + 1 for nested domains
682	!-- Instead it is set to the values of nzt (see routine vnest_boundary_conds
683	!-- or pmci_interp_tril_t) BEFORE calling the pressure solver. To avoid jumps
684	!-- while plotting profiles w at the top has to be set to the values in the
685	!-- height nzt after above modifications. Hint: w level nzt+1 does not impact
686	!-- results.
687	IF ( child_domain .OR. coupling_mode == 'vnested_fine' ) THEN
688	w(nzt+1,:,:) = w(nzt,:,:)
689	ENDIF
690
691	!
692	!-- Sum up the volume flow through the right and north boundary
693	IF ( conserve_volume_flow .AND. bc_lr_cyc .AND. bc_ns_cyc .AND. &
694	nxr == nx ) THEN
695
696	!$OMP PARALLEL PRIVATE (j,k)
697	!$OMP DO
698	DO j = nys, nyn
699	!$OMP CRITICAL
700	DO k = nzb+1, nzt
701	volume_flow_l(1) = volume_flow_l(1) + u(k,j,nxr) * dzw(k) &
702	* MERGE( 1.0_wp, 0.0_wp, &
703	BTEST( wall_flags_total_0(k,j,nxr), 1 )&
704	)
705	ENDDO
706	!$OMP END CRITICAL
707	ENDDO
708	!$OMP END PARALLEL
709
710	ENDIF
711
712	IF ( conserve_volume_flow .AND. bc_ns_cyc .AND. bc_lr_cyc .AND. &
713	nyn == ny ) THEN
714
715	!$OMP PARALLEL PRIVATE (i,k)
716	!$OMP DO
717	DO i = nxl, nxr
718	!$OMP CRITICAL
719	DO k = nzb+1, nzt
720	volume_flow_l(2) = volume_flow_l(2) + v(k,nyn,i) * dzw(k) &
721	* MERGE( 1.0_wp, 0.0_wp, &
722	BTEST( wall_flags_total_0(k,nyn,i), 2 )&
723	)
724	ENDDO
725	!$OMP END CRITICAL
726	ENDDO
727	!$OMP END PARALLEL
728
729	ENDIF
730
731	!
732	!-- Conserve the volume flow
733	IF ( conserve_volume_flow .AND. ( bc_lr_cyc .AND. bc_ns_cyc ) ) THEN
734
735	#if defined( __parallel )
736	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
737	CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 2, MPI_REAL, &
738	MPI_SUM, comm2d, ierr )
739	#else
740	volume_flow = volume_flow_l
741	#endif
742
743	volume_flow_offset(1:2) = ( volume_flow_initial(1:2) - volume_flow(1:2) ) / &
744	volume_flow_area(1:2)
745
746	!$OMP PARALLEL PRIVATE (i,j,k)
747	!$OMP DO
748	DO i = nxl, nxr
749	DO j = nys, nyn
750	DO k = nzb+1, nzt
751	u(k,j,i) = u(k,j,i) + volume_flow_offset(1) &
752	* MERGE( 1.0_wp, 0.0_wp, &
753	BTEST( wall_flags_total_0(k,j,i), 1 ) &
754	)
755	ENDDO
756	DO k = nzb+1, nzt
757	v(k,j,i) = v(k,j,i) + volume_flow_offset(2) &
758	* MERGE( 1.0_wp, 0.0_wp, &
759	BTEST( wall_flags_total_0(k,j,i), 2 ) &
760	)
761	ENDDO
762	ENDDO
763	ENDDO
764
765	!$OMP END PARALLEL
766
767	ENDIF
768
769	!
770	!-- Exchange of boundaries for the velocities
771	CALL exchange_horiz( u, nbgp )
772	CALL exchange_horiz( v, nbgp )
773	CALL exchange_horiz( w, nbgp )
774
775	!
776	!-- Compute the divergence of the corrected velocity field,
777	!-- A possible PE-sum is computed in flow_statistics. Carry out computation
778	!-- only at last Runge-Kutta step.
779	IF ( intermediate_timestep_count == intermediate_timestep_count_max .OR. &
780	intermediate_timestep_count == 0 ) THEN
781	CALL cpu_log( log_point_s(1), 'divergence', 'start' )
782	sums_divnew_l = 0.0_wp
783
784	!
785	!-- d must be reset to zero because it can contain nonzero values below the
786	!-- topography
787	IF ( topography /= 'flat' ) d = 0.0_wp
788
789	localsum = 0.0_wp
790	threadsum = 0.0_wp
791
792	!$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
793	#if defined( __ibm )
794	!$OMP DO SCHEDULE( STATIC )
795	DO i = nxl, nxr
796	DO j = nys, nyn
797	DO k = nzb+1, nzt
798	d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx + &
799	( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy + &
800	( w(k,j,i) * rho_air_zw(k) - &
801	w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k) &
802	) * MERGE( 1.0_wp, 0.0_wp, &
803	BTEST( wall_flags_total_0(k,j,i), 0 ) &
804	)
805	ENDDO
806	DO k = nzb+1, nzt
807	threadsum = threadsum + ABS( d(k,j,i) )
808	ENDDO
809	ENDDO
810	ENDDO
811	#else
812	!$OMP DO SCHEDULE( STATIC )
813	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
814	!$ACC PRESENT(u, v, w, rho_air, rho_air_zw, ddzw, wall_flags_total_0) &
815	!$ACC PRESENT(d)
816	DO i = nxl, nxr
817	DO j = nys, nyn
818	DO k = nzb+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	) * MERGE( 1.0_wp, 0.0_wp, &
824	BTEST( wall_flags_total_0(k,j,i), 0 ) &
825	)
826	ENDDO
827	ENDDO
828	ENDDO
829	!
830	!-- Compute possible PE-sum of divergences for flow_statistics
831	!$OMP DO SCHEDULE( STATIC )
832	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
833	!$ACC REDUCTION(+:threadsum) COPY(threadsum) &
834	!$ACC PRESENT(d)
835	DO i = nxl, nxr
836	DO j = nys, nyn
837	DO k = nzb+1, nzt
838	threadsum = threadsum + ABS( d(k,j,i) )
839	ENDDO
840	ENDDO
841	ENDDO
842	#endif
843
844	localsum = localsum + threadsum
845	!$OMP END PARALLEL
846
847	!
848	!-- For completeness, set the divergence sum of all statistic regions to those
849	!-- of the total domain
850	sums_divnew_l(0:statistic_regions) = localsum
851
852	CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
853
854	ENDIF
855
856	CALL cpu_log( log_point(8), 'pres', 'stop' )
857
858	END SUBROUTINE pres

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |