Home

Context Navigation

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

Last change on this file since 1802 was 1800, checked in by gronemeier, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 24.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |