Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |