Home

Context Navigation

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

Last change on this file since 3318 was 3183, checked in by suehring, 6 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 31.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |