Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |