Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |