Home

Context Navigation

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

Last change on this file since 1575 was 1575, checked in by raasch, 10 years ago
optimized multigrid method installed, new parameter seed_follows_topography for particle release, small adjustment in subjob for HLRN
Property svn:keywords set to `Id`
File size: 24.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |