Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |