Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |