Home

Context Navigation

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

Last change on this file since 2235 was 2233, checked in by suehring, 7 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 30.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |