Home

Context Navigation

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

Last change on this file since 2039 was 2038, checked in by knoop, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 27.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |