Home

Context Navigation

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

Last change on this file since 4649 was 4649, checked in by raasch, 4 years ago
files re-formatted to follow the PALM coding standard
Property svn:keywords set to `Id`
File size: 33.4 KB

Rev	Line
[1682]	1	!> @file pres.f90
[4649]	2	!--------------------------------------------------------------------------------------------------!
[2696]	3	! This file is part of the PALM model system.
[1036]	4	!
[4649]	5	! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
	6	! Public License as published by the Free Software Foundation, either version 3 of the License, or
	7	! (at your option) any later version.
[1036]	8	!
[4649]	9	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
	10	! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
	11	! Public License for more details.
[1036]	12	!
[4649]	13	! You should have received a copy of the GNU General Public License along with PALM. If not, see
	14	! <http://www.gnu.org/licenses/>.
[1036]	15	!
[4360]	16	! Copyright 1997-2020 Leibniz Universitaet Hannover
[4649]	17	!--------------------------------------------------------------------------------------------------!
[1036]	18	!
[4649]	19	!
[484]	20	! Current revisions:
[4649]	21	! -----------------
	22	!
	23	!
[1919]	24	! Former revisions:
	25	! -----------------
[3016]	26	! $Id: pres.f90 4649 2020-08-25 12:11:17Z raasch $
[4649]	27	! File re-formatted to follow the PALM coding standard
	28	!
	29	!
	30	! 4457 2020-03-11 14:20:43Z raasch
[4457]	31	! use statement for exchange horiz added
[4649]	32	!
[4457]	33	! 4360 2020-01-07 11:25:50Z suehring
[4346]	34	! Introduction of wall_flags_total_0, which currently sets bits based on static
	35	! topography information used in wall_flags_static_0
[4649]	36	!
[4346]	37	! 4329 2019-12-10 15:46:36Z motisi
[4329]	38	! Renamed wall_flags_0 to wall_flags_static_0
[4649]	39	!
[4329]	40	! 4182 2019-08-22 15:20:23Z scharf
[4182]	41	! Corrected "Former revisions" section
[4649]	42	!
[4182]	43	! 4015 2019-06-05 13:25:35Z raasch
[4015]	44	! variable child_domain_nvn eliminated
[4649]	45	!
[4015]	46	! 3849 2019-04-01 16:35:16Z knoop
[3634]	47	! OpenACC port for SPEC
[1919]	48	!
[4182]	49	! Revision 1.1 1997/07/24 11:24:44 raasch
	50	! Initial revision
	51	!
	52	!
[4649]	53	!--------------------------------------------------------------------------------------------------!
[1]	54	! Description:
	55	! ------------
[4649]	56	!> Compute the divergence of the provisional velocity field. Solve the Poisson equation for the
	57	!> perturbation pressure. Compute the final velocities using this perturbation pressure. Compute the
	58	!> remaining divergence.
	59	!--------------------------------------------------------------------------------------------------!
[1682]	60	SUBROUTINE pres
[1]	61
[1320]	62
[4649]	63	USE arrays_3d, &
	64	ONLY: d, &
	65	ddzu, &
	66	ddzu_pres, &
	67	ddzw, &
	68	dzw, &
	69	p, &
	70	p_loc, &
	71	rho_air, &
	72	rho_air_zw, &
	73	tend, &
	74	u, &
	75	v, &
	76	w
[1320]	77
[4649]	78	USE control_parameters, &
	79	ONLY: bc_lr_cyc, &
	80	bc_ns_cyc, &
	81	bc_radiation_l, &
	82	bc_radiation_n, &
	83	bc_radiation_r, &
	84	bc_radiation_s, &
	85	child_domain, &
	86	conserve_volume_flow, &
	87	coupling_mode, &
	88	dt_3d, &
	89	gathered_size, &
	90	ibc_p_b, &
	91	ibc_p_t, &
	92	intermediate_timestep_count, &
	93	intermediate_timestep_count_max, &
	94	mg_switch_to_pe0_level, &
	95	nesting_offline, &
	96	psolver, &
	97	subdomain_size, &
	98	topography, &
	99	volume_flow, &
	100	volume_flow_area, &
	101	volume_flow_initial
[1320]	102
[4649]	103	USE cpulog, &
	104	ONLY: cpu_log, &
	105	log_point, &
	106	log_point_s
	107
	108	USE exchange_horiz_mod, &
[4457]	109	ONLY: exchange_horiz
	110
[4649]	111	USE grid_variables, &
	112	ONLY: ddx, &
	113	ddy
[1320]	114
[4649]	115	USE indices, &
	116	ONLY: nbgp, &
	117	ngp_2dh_outer, &
	118	nx, &
	119	nxl, &
	120	nxlg, &
	121	nxl_mg, &
	122	nxr, &
	123	nxrg, &
	124	nxr_mg, &
	125	ny, &
	126	nys, &
	127	nysg, &
	128	nys_mg, &
	129	nyn, &
	130	nyng, &
	131	nyn_mg, &
	132	nzb, &
	133	nzt, &
	134	nzt_mg, &
[4346]	135	wall_flags_total_0
[1320]	136
	137	USE kinds
	138
[1]	139	USE pegrid
	140
[4649]	141	USE pmc_interface, &
	142	ONLY: nesting_mode
	143
	144	USE poisfft_mod, &
[1320]	145	ONLY: poisfft
	146
[1575]	147	USE poismg_mod
	148
[2696]	149	USE poismg_noopt_mod
	150
[4649]	151	USE statistics, &
	152	ONLY: statistic_regions, &
	153	sums_divnew_l, &
	154	sums_divold_l, &
	155	weight_pres, &
[1320]	156	weight_substep
	157
[4649]	158	USE surface_mod, &
[2232]	159	ONLY : bc_h
	160
[1]	161	IMPLICIT NONE
	162
[4649]	163	INTEGER(iwp) :: i !<
	164	INTEGER(iwp) :: j !<
	165	INTEGER(iwp) :: k !<
	166	INTEGER(iwp) :: m !<
[1]	167
[4649]	168	REAL(wp) :: ddt_3d !<
	169	REAL(wp) :: d_weight_pres !<
	170	REAL(wp) :: localsum !<
	171	REAL(wp) :: threadsum !<
	172	REAL(wp) :: weight_pres_l !<
	173	REAL(wp) :: weight_substep_l !<
[1]	174
[1762]	175	REAL(wp), DIMENSION(1:3) :: volume_flow_l !<
	176	REAL(wp), DIMENSION(1:3) :: volume_flow_offset !<
[1682]	177	REAL(wp), DIMENSION(1:nzt) :: w_l !<
	178	REAL(wp), DIMENSION(1:nzt) :: w_l_l !<
[1]	179
	180
	181	CALL cpu_log( log_point(8), 'pres', 'start' )
	182
[1918]	183	!
	184	!-- Calculate quantities to be used locally
[1342]	185	ddt_3d = 1.0_wp / dt_3d
[1918]	186	IF ( intermediate_timestep_count == 0 ) THEN
	187	!
	188	!-- If pres is called before initial time step
[1929]	189	weight_pres_l = 1.0_wp
	190	d_weight_pres = 1.0_wp
	191	weight_substep_l = 1.0_wp
[1918]	192	ELSE
[1929]	193	weight_pres_l = weight_pres(intermediate_timestep_count)
	194	d_weight_pres = 1.0_wp / weight_pres(intermediate_timestep_count)
	195	weight_substep_l = weight_substep(intermediate_timestep_count)
[1918]	196	ENDIF
[85]	197
[1]	198	!
[707]	199	!-- Multigrid method expects array d to have one ghost layer.
[4649]	200	!--
[1575]	201	IF ( psolver(1:9) == 'multigrid' ) THEN
[4649]	202
[1]	203	DEALLOCATE( d )
[4649]	204	ALLOCATE( d(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
[707]	205
	206	!
[4649]	207	!-- Since p is later used to hold the weighted average of the substeps, it cannot be used in the
	208	!-- iterative solver. Therefore, its initial value is stored on p_loc, which is then iteratively
	209	!-- advanced in every substep.
[1762]	210	IF ( intermediate_timestep_count <= 1 ) THEN
[707]	211	DO i = nxl-1, nxr+1
	212	DO j = nys-1, nyn+1
	213	DO k = nzb, nzt+1
	214	p_loc(k,j,i) = p(k,j,i)
	215	ENDDO
	216	ENDDO
	217	ENDDO
[667]	218	ENDIF
[4649]	219
[1918]	220	ELSEIF ( psolver == 'sor' .AND. intermediate_timestep_count <= 1 ) THEN
[707]	221
	222	!
[4649]	223	!-- Since p is later used to hold the weighted average of the substeps, it cannot be used in the
	224	!-- iterative solver. Therefore, its initial value is stored on p_loc, which is then iteratively
	225	!-- advanced in every substep.
[707]	226	p_loc = p
	227
[1]	228	ENDIF
	229
	230	!
[4649]	231	!-- Conserve the volume flow at the outflow in case of non-cyclic lateral boundary conditions
	232	!-- WARNING: so far, this conservation does not work at the left/south boundary if the topography at
	233	!-- the inflow differs from that at the outflow! For this case, volume_flow_area needs
	234	!-- adjustment!
[106]	235	!
	236	!-- Left/right
[4649]	237	IF ( conserve_volume_flow .AND. ( bc_radiation_l .OR. bc_radiation_r ) ) THEN
[680]	238
[1342]	239	volume_flow(1) = 0.0_wp
	240	volume_flow_l(1) = 0.0_wp
[106]	241
[3182]	242	IF ( bc_radiation_l ) THEN
[106]	243	i = 0
[3182]	244	ELSEIF ( bc_radiation_r ) THEN
[106]	245	i = nx+1
	246	ENDIF
	247
	248	DO j = nys, nyn
	249	!
	250	!-- Sum up the volume flow through the south/north boundary
[2232]	251	DO k = nzb+1, nzt
[4649]	252	volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k) &
	253	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
[106]	254	ENDDO
	255	ENDDO
	256
[4649]	257	#if defined( __parallel )
[680]	258	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
[4649]	259	CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 1, MPI_REAL, MPI_SUM, comm1dy, ierr )
[106]	260	#else
[4649]	261	volume_flow = volume_flow_l
[106]	262	#endif
[4649]	263	volume_flow_offset(1) = ( volume_flow_initial(1) - volume_flow(1) ) / volume_flow_area(1)
[106]	264
[667]	265	DO j = nysg, nyng
[2232]	266	DO k = nzb+1, nzt
[4649]	267	u(k,j,i) = u(k,j,i) + volume_flow_offset(1) &
	268	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
[106]	269	ENDDO
	270	ENDDO
	271
	272	ENDIF
	273
	274	!
	275	!-- South/north
[3182]	276	IF ( conserve_volume_flow .AND. ( bc_radiation_n .OR. bc_radiation_s ) ) THEN
[106]	277
[1342]	278	volume_flow(2) = 0.0_wp
	279	volume_flow_l(2) = 0.0_wp
[75]	280
[3182]	281	IF ( bc_radiation_s ) THEN
[106]	282	j = 0
[3182]	283	ELSEIF ( bc_radiation_n ) THEN
[75]	284	j = ny+1
[106]	285	ENDIF
	286
	287	DO i = nxl, nxr
[75]	288	!
[106]	289	!-- Sum up the volume flow through the south/north boundary
[2232]	290	DO k = nzb+1, nzt
[4649]	291	volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k) &
	292	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
[75]	293	ENDDO
[106]	294	ENDDO
	295
[4649]	296	#if defined( __parallel )
[680]	297	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
[4649]	298	CALL MPI_ALLREDUCE( volume_flow_l(2), volume_flow(2), 1, MPI_REAL, MPI_SUM, comm1dx, ierr )
[75]	299	#else
[4649]	300	volume_flow = volume_flow_l
[75]	301	#endif
[4649]	302	volume_flow_offset(2) = ( volume_flow_initial(2) - volume_flow(2) ) / volume_flow_area(2)
[75]	303
[667]	304	DO i = nxlg, nxrg
[2232]	305	DO k = nzb+1, nzt
[4649]	306	v(k,j,i) = v(k,j,i) + volume_flow_offset(2) &
	307	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
[75]	308	ENDDO
[106]	309	ENDDO
[75]	310
	311	ENDIF
	312
[76]	313	!
[4015]	314	!-- Remove mean vertical velocity in case that Neumann conditions are used both at bottom and top
[4649]	315	!-- boundary. With Neumann conditions at both vertical boundaries, the solver cannot remove mean
	316	!-- vertical velocities. They should be removed, because incompressibility requires that the
	317	!-- vertical gradient of vertical velocity is zero. Since w=0 at the solid surface, it must be zero
	318	!-- everywhere.
[4015]	319	!-- This must not be done in case of a 3d-nesting child domain, because a mean vertical velocity
	320	!-- can physically exist in such a domain.
	321	!-- Also in case of offline nesting, mean vertical velocities may exist (and must not be removed),
	322	!-- caused by horizontal divergence/convergence of the large scale flow that is prescribed at the
	323	!-- side boundaries.
[4649]	324	!-- The removal cannot be done before the first initial time step because ngp_2dh_outer is not yet
	325	!-- known then.
[4015]	326	IF ( ibc_p_b == 1 .AND. ibc_p_t == 1 .AND. .NOT. nesting_offline &
	327	.AND. .NOT. ( child_domain .AND. nesting_mode /= 'vertical' ) &
[4649]	328	.AND. intermediate_timestep_count /= 0 ) THEN
[1918]	329	w_l = 0.0_wp; w_l_l = 0.0_wp
	330	DO i = nxl, nxr
	331	DO j = nys, nyn
[2232]	332	DO k = nzb+1, nzt
[4015]	333	w_l_l(k) = w_l_l(k) + w(k,j,i) &
[4649]	334	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
[76]	335	ENDDO
	336	ENDDO
[1918]	337	ENDDO
[4649]	338	#if defined( __parallel )
[1918]	339	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[4015]	340	CALL MPI_ALLREDUCE( w_l_l(1), w_l(1), nzt, MPI_REAL, MPI_SUM, comm2d, ierr )
[76]	341	#else
[1918]	342	w_l = w_l_l
[76]	343	#endif
[1918]	344	DO k = 1, nzt
	345	w_l(k) = w_l(k) / ngp_2dh_outer(k,0)
	346	ENDDO
	347	DO i = nxlg, nxrg
	348	DO j = nysg, nyng
[2232]	349	DO k = nzb+1, nzt
[4015]	350	w(k,j,i) = w(k,j,i) - w_l(k) &
[4649]	351	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
[76]	352	ENDDO
	353	ENDDO
[1918]	354	ENDDO
[76]	355	ENDIF
[75]	356
	357	!
[1]	358	!-- Compute the divergence of the provisional velocity field.
	359	CALL cpu_log( log_point_s(1), 'divergence', 'start' )
	360
[1575]	361	IF ( psolver(1:9) == 'multigrid' ) THEN
[2232]	362	!$OMP PARALLEL DO SCHEDULE( STATIC ) PRIVATE (i,j,k)
[1]	363	DO i = nxl-1, nxr+1
	364	DO j = nys-1, nyn+1
	365	DO k = nzb, nzt+1
[1342]	366	d(k,j,i) = 0.0_wp
[1]	367	ENDDO
	368	ENDDO
	369	ENDDO
	370	ELSE
[2232]	371	!$OMP PARALLEL DO SCHEDULE( STATIC ) PRIVATE (i,j,k)
[3634]	372	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
	373	!$ACC PRESENT(d)
[1003]	374	DO i = nxl, nxr
	375	DO j = nys, nyn
	376	DO k = nzb+1, nzt
[1342]	377	d(k,j,i) = 0.0_wp
[1]	378	ENDDO
	379	ENDDO
	380	ENDDO
	381	ENDIF
	382
[1342]	383	localsum = 0.0_wp
	384	threadsum = 0.0_wp
[1]	385
	386	#if defined( __ibm )
	387	!$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
	388	!$OMP DO SCHEDULE( STATIC )
	389	DO i = nxl, nxr
	390	DO j = nys, nyn
[2232]	391	DO k = nzb+1, nzt
[4649]	392	d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx + &
	393	( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy + &
	394	( w(k,j,i) * rho_air_zw(k) - w(k-1,j,i) * rho_air_zw(k-1) ) &
	395	* ddzw(k) ) * ddt_3d * d_weight_pres &
	396	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
[1]	397	ENDDO
	398	!
	399	!-- Compute possible PE-sum of divergences for flow_statistics
[2232]	400	DO k = nzb+1, nzt
[4649]	401	threadsum = threadsum + ABS( d(k,j,i) ) &
	402	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
[1]	403	ENDDO
	404
	405	ENDDO
	406	ENDDO
	407
[4649]	408	IF ( intermediate_timestep_count == intermediate_timestep_count_max .OR. &
[1918]	409	intermediate_timestep_count == 0 ) THEN
	410	localsum = localsum + threadsum * dt_3d * weight_pres_l
[1908]	411	ENDIF
[1]	412	!$OMP END PARALLEL
	413	#else
	414
[1111]	415	!$OMP PARALLEL PRIVATE (i,j,k)
	416	!$OMP DO SCHEDULE( STATIC )
[3634]	417	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
[4346]	418	!$ACC PRESENT(u, v, w, rho_air, rho_air_zw, ddzw, wall_flags_total_0) &
[3634]	419	!$ACC PRESENT(d)
[1111]	420	DO i = nxl, nxr
	421	DO j = nys, nyn
	422	DO k = 1, nzt
[4649]	423	d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx + &
	424	( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy + &
	425	( w(k,j,i) * rho_air_zw(k) - w(k-1,j,i) * rho_air_zw(k-1) ) &
	426	* ddzw(k) ) * ddt_3d * d_weight_pres &
	427	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
[1]	428	ENDDO
	429	ENDDO
[1111]	430	ENDDO
	431	!$OMP END PARALLEL
[1]	432
	433	!
[4649]	434	!-- Compute possible PE-sum of divergences for flow_statistics. Carry out computation only at last
	435	!-- Runge-Kutta substep.
	436	IF ( intermediate_timestep_count == intermediate_timestep_count_max .OR. &
[1918]	437	intermediate_timestep_count == 0 ) THEN
[1908]	438	!$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
	439	!$OMP DO SCHEDULE( STATIC )
[3634]	440	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
	441	!$ACC REDUCTION(+:threadsum) COPY(threadsum) &
	442	!$ACC PRESENT(d)
[1908]	443	DO i = nxl, nxr
	444	DO j = nys, nyn
	445	DO k = nzb+1, nzt
	446	threadsum = threadsum + ABS( d(k,j,i) )
	447	ENDDO
[1]	448	ENDDO
	449	ENDDO
[1918]	450	localsum = localsum + threadsum * dt_3d * weight_pres_l
[1908]	451	!$OMP END PARALLEL
	452	ENDIF
[1]	453	#endif
	454
	455	!
[4649]	456	!-- For completeness, set the divergence sum of all statistic regions to those of the total domain
	457	IF ( intermediate_timestep_count == intermediate_timestep_count_max .OR. &
[1918]	458	intermediate_timestep_count == 0 ) THEN
[1908]	459	sums_divold_l(0:statistic_regions) = localsum
[1918]	460	ENDIF
[1]	461
	462	CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
	463
	464	!
	465	!-- Compute the pressure perturbation solving the Poisson equation
[1575]	466	IF ( psolver == 'poisfft' ) THEN
[1]	467
	468	!
	469	!-- Solve Poisson equation via FFT and solution of tridiagonal matrices
[1575]	470	CALL poisfft( d )
[1212]	471
[1]	472	!
[4649]	473	!-- Store computed perturbation pressure and set boundary condition in z-direction
[2232]	474	!$OMP PARALLEL DO PRIVATE (i,j,k)
[3634]	475	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
	476	!$ACC PRESENT(d, tend)
[1]	477	DO i = nxl, nxr
	478	DO j = nys, nyn
	479	DO k = nzb+1, nzt
	480	tend(k,j,i) = d(k,j,i)
	481	ENDDO
	482	ENDDO
	483	ENDDO
	484
	485	!
	486	!-- Bottom boundary:
[4649]	487	!-- This condition is only required for internal output. The pressure gradient
	488	!-- (dp(nzb+1)-dp(nzb))/dz is not used anywhere else.
[1]	489	IF ( ibc_p_b == 1 ) THEN
	490	!
[4649]	491	!-- Neumann (dp/dz = 0). Using surface data type, first for non-natural surfaces, then for
	492	!-- natural and urban surfaces
[2232]	493	!-- Upward facing
	494	!$OMP PARALLEL DO PRIVATE( i, j, k )
[3634]	495	!$ACC PARALLEL LOOP PRIVATE(i, j, k) &
	496	!$ACC PRESENT(bc_h, tend)
[2232]	497	DO m = 1, bc_h(0)%ns
[2298]	498	i = bc_h(0)%i(m)
	499	j = bc_h(0)%j(m)
	500	k = bc_h(0)%k(m)
[2232]	501	tend(k-1,j,i) = tend(k,j,i)
[1]	502	ENDDO
[2232]	503	!
	504	!-- Downward facing
	505	!$OMP PARALLEL DO PRIVATE( i, j, k )
[3634]	506	!$ACC PARALLEL LOOP PRIVATE(i, j, k) &
	507	!$ACC PRESENT(bc_h, tend)
[2232]	508	DO m = 1, bc_h(1)%ns
[2298]	509	i = bc_h(1)%i(m)
	510	j = bc_h(1)%j(m)
	511	k = bc_h(1)%k(m)
[2232]	512	tend(k+1,j,i) = tend(k,j,i)
	513	ENDDO
[1]	514
	515	ELSE
	516	!
[4649]	517	!-- Dirichlet. Using surface data type, first for non-natural surfaces, then for natural and
	518	!-- urban surfaces
[2232]	519	!-- Upward facing
	520	!$OMP PARALLEL DO PRIVATE( i, j, k )
	521	DO m = 1, bc_h(0)%ns
[2298]	522	i = bc_h(0)%i(m)
	523	j = bc_h(0)%j(m)
	524	k = bc_h(0)%k(m)
[2232]	525	tend(k-1,j,i) = 0.0_wp
[1]	526	ENDDO
[2232]	527	!
	528	!-- Downward facing
	529	!$OMP PARALLEL DO PRIVATE( i, j, k )
	530	DO m = 1, bc_h(1)%ns
[2298]	531	i = bc_h(1)%i(m)
	532	j = bc_h(1)%j(m)
	533	k = bc_h(1)%k(m)
[2232]	534	tend(k+1,j,i) = 0.0_wp
	535	ENDDO
[1]	536
	537	ENDIF
	538
	539	!
	540	!-- Top boundary
	541	IF ( ibc_p_t == 1 ) THEN
	542	!
	543	!-- Neumann
[2232]	544	!$OMP PARALLEL DO PRIVATE (i,j,k)
[667]	545	DO i = nxlg, nxrg
	546	DO j = nysg, nyng
[1]	547	tend(nzt+1,j,i) = tend(nzt,j,i)
	548	ENDDO
	549	ENDDO
	550
	551	ELSE
	552	!
	553	!-- Dirichlet
[2232]	554	!$OMP PARALLEL DO PRIVATE (i,j,k)
[3634]	555	!$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
	556	!$ACC PRESENT(tend)
[667]	557	DO i = nxlg, nxrg
	558	DO j = nysg, nyng
[1342]	559	tend(nzt+1,j,i) = 0.0_wp
[1]	560	ENDDO
	561	ENDDO
	562
	563	ENDIF
	564
	565	!
	566	!-- Exchange boundaries for p
[667]	567	CALL exchange_horiz( tend, nbgp )
[4649]	568
[1]	569	ELSEIF ( psolver == 'sor' ) THEN
	570
	571	!
[4649]	572	!-- Solve Poisson equation for perturbation pressure using SOR-Red/Black scheme
[707]	573	CALL sor( d, ddzu_pres, ddzw, p_loc )
	574	tend = p_loc
[1]	575
[1575]	576	ELSEIF ( psolver(1:9) == 'multigrid' ) THEN
[1]	577
	578	!
[4649]	579	!-- Solve Poisson equation for perturbation pressure using Multigrid scheme, array tend is used
	580	!-- to store the residuals.
[778]	581
[4649]	582	!-- If the number of grid points of the gathered grid, which is collected on PE0, is larger than
	583	!-- the number of grid points of an PE, than array tend will be enlarged.
[778]	584	IF ( gathered_size > subdomain_size ) THEN
	585	DEALLOCATE( tend )
[4649]	586	ALLOCATE( tend(nzb:nzt_mg(mg_switch_to_pe0_level)+1,nys_mg( &
	587	mg_switch_to_pe0_level)-1:nyn_mg(mg_switch_to_pe0_level)+1, &
	588	nxl_mg(mg_switch_to_pe0_level)-1:nxr_mg(mg_switch_to_pe0_level)+1) )
[778]	589	ENDIF
	590
[1575]	591	IF ( psolver == 'multigrid' ) THEN
	592	CALL poismg( tend )
	593	ELSE
[1931]	594	CALL poismg_noopt( tend )
[1575]	595	ENDIF
[707]	596
[778]	597	IF ( gathered_size > subdomain_size ) THEN
	598	DEALLOCATE( tend )
	599	ALLOCATE( tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	600	ENDIF
	601
[1]	602	!
[4649]	603	!-- Restore perturbation pressure on tend because this array is used further below to correct the
	604	!-- velocity fields
[707]	605	DO i = nxl-1, nxr+1
	606	DO j = nys-1, nyn+1
	607	DO k = nzb, nzt+1
	608	tend(k,j,i) = p_loc(k,j,i)
	609	ENDDO
	610	ENDDO
	611	ENDDO
[667]	612
[1]	613	ENDIF
	614
	615	!
[707]	616	!-- Store perturbation pressure on array p, used for pressure data output.
	617	!-- Ghost layers are added in the output routines (except sor-method: see below)
[1918]	618	IF ( intermediate_timestep_count <= 1 ) THEN
[707]	619	!$OMP PARALLEL PRIVATE (i,j,k)
	620	!$OMP DO
[3634]	621	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
	622	!$ACC PRESENT(p, tend)
[707]	623	DO i = nxl-1, nxr+1
	624	DO j = nys-1, nyn+1
	625	DO k = nzb, nzt+1
[4649]	626	p(k,j,i) = tend(k,j,i) * weight_substep_l
[673]	627	ENDDO
[1]	628	ENDDO
[707]	629	ENDDO
	630	!$OMP END PARALLEL
	631
[1762]	632	ELSEIF ( intermediate_timestep_count > 1 ) THEN
[707]	633	!$OMP PARALLEL PRIVATE (i,j,k)
	634	!$OMP DO
[3634]	635	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
	636	!$ACC PRESENT(p, tend)
[707]	637	DO i = nxl-1, nxr+1
	638	DO j = nys-1, nyn+1
	639	DO k = nzb, nzt+1
[4649]	640	p(k,j,i) = p(k,j,i) + tend(k,j,i) * weight_substep_l
[673]	641	ENDDO
	642	ENDDO
[707]	643	ENDDO
	644	!$OMP END PARALLEL
	645
	646	ENDIF
[4649]	647
[707]	648	!
	649	!-- SOR-method needs ghost layers for the next timestep
	650	IF ( psolver == 'sor' ) CALL exchange_horiz( p, nbgp )
[682]	651
[1]	652	!
[4649]	653	!-- Correction of the provisional velocities with the current perturbation pressure just computed
[709]	654	IF ( conserve_volume_flow .AND. ( bc_lr_cyc .OR. bc_ns_cyc ) ) THEN
[1342]	655	volume_flow_l(1) = 0.0_wp
	656	volume_flow_l(2) = 0.0_wp
[1]	657	ENDIF
[3347]	658	!
[4649]	659	!-- Add pressure gradients to the velocity components. Note, the loops are running over the entire
	660	!-- model domain, even though, in case of non-cyclic boundaries u- and v-component are not
	661	!-- prognostic at i=0 and j=0, respectiveley. However, in case of Dirichlet boundary conditions for
	662	!-- the velocities, zero-gradient conditions for the pressure are set, so that no modification is
	663	!-- imposed at the boundaries.
[1]	664	!$OMP PARALLEL PRIVATE (i,j,k)
	665	!$OMP DO
[3634]	666	!$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k) &
[4346]	667	!$ACC PRESENT(u, v, w, tend, ddzu, wall_flags_total_0)
[4649]	668	DO i = nxl, nxr
[1]	669	DO j = nys, nyn
[2118]	670
[2232]	671	DO k = nzb+1, nzt
[4649]	672	w(k,j,i) = w(k,j,i) - dt_3d * ( tend(k+1,j,i) - tend(k,j,i) ) * ddzu(k+1) &
	673	* weight_pres_l &
	674	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
[1]	675	ENDDO
[2118]	676
[2232]	677	DO k = nzb+1, nzt
[4649]	678	u(k,j,i) = u(k,j,i) - dt_3d * ( tend(k,j,i) - tend(k,j,i-1) ) * ddx * weight_pres_l &
	679	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
[1]	680	ENDDO
[2118]	681
[2232]	682	DO k = nzb+1, nzt
[4649]	683	v(k,j,i) = v(k,j,i) - dt_3d * ( tend(k,j,i) - tend(k,j-1,i) ) * ddy * weight_pres_l &
	684	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
	685	ENDDO
[1]	686
	687	ENDDO
	688	ENDDO
	689	!$OMP END PARALLEL
[1113]	690
	691	!
[4649]	692	!-- The vertical velocity is not set to zero at nzt + 1 for nested domains. Instead it is set to the
	693	!-- values of nzt (see routine vnest_boundary_conds or pmci_interp_tril_t) BEFORE calling the
	694	!-- pressure solver. To avoid jumps while plotting profiles, w at the top has to be set to the
	695	!-- values in height nzt after above modifications. Hint: w level nzt+1 does not impact results.
	696	IF ( child_domain .OR. coupling_mode == 'vnested_fine' ) THEN
[3012]	697	w(nzt+1,:,:) = w(nzt,:,:)
	698	ENDIF
	699
	700	!
[1113]	701	!-- Sum up the volume flow through the right and north boundary
[4649]	702	IF ( conserve_volume_flow .AND. bc_lr_cyc .AND. bc_ns_cyc .AND. nxr == nx ) THEN
[1113]	703
	704	!$OMP PARALLEL PRIVATE (j,k)
	705	!$OMP DO
	706	DO j = nys, nyn
	707	!$OMP CRITICAL
[2232]	708	DO k = nzb+1, nzt
[4649]	709	volume_flow_l(1) = volume_flow_l(1) + u(k,j,nxr) * dzw(k) &
	710	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,nxr), 1 ) )
[1113]	711	ENDDO
	712	!$OMP END CRITICAL
	713	ENDDO
	714	!$OMP END PARALLEL
	715
	716	ENDIF
	717
[4649]	718	IF ( conserve_volume_flow .AND. bc_ns_cyc .AND. bc_lr_cyc .AND. nyn == ny ) THEN
[1113]	719
	720	!$OMP PARALLEL PRIVATE (i,k)
	721	!$OMP DO
	722	DO i = nxl, nxr
	723	!$OMP CRITICAL
[2232]	724	DO k = nzb+1, nzt
[4649]	725	volume_flow_l(2) = volume_flow_l(2) + v(k,nyn,i) * dzw(k) &
	726	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,nyn,i), 2 ) )
[1113]	727	ENDDO
	728	!$OMP END CRITICAL
	729	ENDDO
	730	!$OMP END PARALLEL
	731
	732	ENDIF
[4649]	733
[1]	734	!
	735	!-- Conserve the volume flow
[707]	736	IF ( conserve_volume_flow .AND. ( bc_lr_cyc .AND. bc_ns_cyc ) ) THEN
[1]	737
[4649]	738	#if defined( __parallel )
[622]	739	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[4649]	740	CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 2, MPI_REAL, MPI_SUM, comm2d, ierr )
[1]	741	#else
[4649]	742	volume_flow = volume_flow_l
	743	#endif
[1]	744
[4649]	745	volume_flow_offset(1:2) = ( volume_flow_initial(1:2) - volume_flow(1:2) ) &
	746	/ volume_flow_area(1:2)
[1]	747
	748	!$OMP PARALLEL PRIVATE (i,j,k)
	749	!$OMP DO
	750	DO i = nxl, nxr
	751	DO j = nys, nyn
[2232]	752	DO k = nzb+1, nzt
[4649]	753	u(k,j,i) = u(k,j,i) + volume_flow_offset(1) &
	754	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
[719]	755	ENDDO
[2232]	756	DO k = nzb+1, nzt
[4649]	757	v(k,j,i) = v(k,j,i) + volume_flow_offset(2) &
	758	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
[667]	759	ENDDO
[1]	760	ENDDO
	761	ENDDO
[667]	762
[1]	763	!$OMP END PARALLEL
	764
	765	ENDIF
	766
	767	!
	768	!-- Exchange of boundaries for the velocities
[667]	769	CALL exchange_horiz( u, nbgp )
	770	CALL exchange_horiz( v, nbgp )
	771	CALL exchange_horiz( w, nbgp )
[1]	772
	773	!
[4649]	774	!-- Compute the divergence of the corrected velocity field.
	775	!-- A possible PE-sum is computed in flow_statistics. Carry out computation only at last
	776	!-- Runge-Kutta step.
	777	IF ( intermediate_timestep_count == intermediate_timestep_count_max .OR. &
[1918]	778	intermediate_timestep_count == 0 ) THEN
[1908]	779	CALL cpu_log( log_point_s(1), 'divergence', 'start' )
	780	sums_divnew_l = 0.0_wp
[1]	781
	782	!
[4649]	783	!-- d must be reset to zero because it can contain nonzero values below the topography
[1908]	784	IF ( topography /= 'flat' ) d = 0.0_wp
[1]	785
[1908]	786	localsum = 0.0_wp
	787	threadsum = 0.0_wp
[1]	788
[1908]	789	!$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
[2073]	790	#if defined( __ibm )
[1908]	791	!$OMP DO SCHEDULE( STATIC )
	792	DO i = nxl, nxr
	793	DO j = nys, nyn
[2232]	794	DO k = nzb+1, nzt
[4649]	795	d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx + &
	796	( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy + &
	797	( w(k,j,i) * rho_air_zw(k) - w(k-1,j,i) * rho_air_zw(k-1) ) &
	798	* ddzw(k) ) &
	799	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
[1908]	800	ENDDO
	801	DO k = nzb+1, nzt
	802	threadsum = threadsum + ABS( d(k,j,i) )
	803	ENDDO
[1]	804	ENDDO
	805	ENDDO
	806	#else
[2073]	807	!$OMP DO SCHEDULE( STATIC )
[3634]	808	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
[4346]	809	!$ACC PRESENT(u, v, w, rho_air, rho_air_zw, ddzw, wall_flags_total_0) &
[3634]	810	!$ACC PRESENT(d)
[1908]	811	DO i = nxl, nxr
	812	DO j = nys, nyn
[2232]	813	DO k = nzb+1, nzt
[2037]	814	d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx + &
	815	( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy + &
[4649]	816	( w(k,j,i) * rho_air_zw(k) - w(k-1,j,i) * rho_air_zw(k-1) ) &
	817	* ddzw(k) ) &
	818	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
[1908]	819	ENDDO
[1113]	820	ENDDO
	821	ENDDO
	822	!
[1908]	823	!-- Compute possible PE-sum of divergences for flow_statistics
[2073]	824	!$OMP DO SCHEDULE( STATIC )
[3634]	825	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
	826	!$ACC REDUCTION(+:threadsum) COPY(threadsum) &
	827	!$ACC PRESENT(d)
[1908]	828	DO i = nxl, nxr
	829	DO j = nys, nyn
	830	DO k = nzb+1, nzt
	831	threadsum = threadsum + ABS( d(k,j,i) )
	832	ENDDO
[1]	833	ENDDO
	834	ENDDO
	835	#endif
[667]	836
[1908]	837	localsum = localsum + threadsum
	838	!$OMP END PARALLEL
[1]	839
	840	!
[4649]	841	!-- For completeness, set the divergence sum of all statistic regions to those of the total
	842	!-- domain
[1908]	843	sums_divnew_l(0:statistic_regions) = localsum
[1]	844
[1908]	845	CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
[1]	846
[1908]	847	ENDIF
	848
[1]	849	CALL cpu_log( log_point(8), 'pres', 'stop' )
	850
[4649]	851	END SUBROUTINE pres

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |