Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |