Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |