Home

Context Navigation

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

Last change on this file since 671 was 668, checked in by suehring, 14 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 20.6 KB

Rev	Line
[1]	1	SUBROUTINE pres
	2
	3	!------------------------------------------------------------------------------!
[484]	4	! Current revisions:
[1]	5	! -----------------
[668]	6	!
	7	! Former revisions:
	8	! -----------------
	9	! $Id: pres.f90 668 2010-12-23 13:22:58Z heinze $
	10	!
	11	! 667 2010-12-23 12:06:00Z suehring/gryschka
[667]	12	! New allocation of tend when ws-scheme and multigrid is used. This is due to
	13	! reasons of perforance of the data_exchange. The same is done with p after
	14	! poismg is called.
	15	! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng when no
	16	! multigrid is used. Calls of exchange_horiz are modified.
	17	! bugfix: After pressure correction no volume flow correction in case of
	18	! non-cyclic boundary conditions
	19	! (has to be done only before pressure correction)
	20	! Call of SOR routine is referenced with ddzu_pres.
	21	!
[623]	22	! 622 2010-12-10 08:08:13Z raasch
	23	! optional barriers included in order to speed up collective operations
	24	!
[198]	25	! 151 2008-03-07 13:42:18Z raasch
	26	! Bugfix in volume flow control for non-cyclic boundary conditions
	27	!
[110]	28	! 106 2007-08-16 14:30:26Z raasch
	29	! Volume flow conservation added for the remaining three outflow boundaries
	30	!
[90]	31	! 85 2007-05-11 09:35:14Z raasch
	32	! Division through dt_3d replaced by multiplication of the inverse.
	33	! For performance optimisation, this is done in the loop calculating the
	34	! divergence instead of using a seperate loop.
	35	!
[77]	36	! 75 2007-03-22 09:54:05Z raasch
[75]	37	! Volume flow control for non-cyclic boundary conditions added (currently only
[76]	38	! for the north boundary!!), 2nd+3rd argument removed from exchange horiz,
	39	! mean vertical velocity is removed in case of Neumann boundary conditions
	40	! both at the bottom and the top
[1]	41	!
[3]	42	! RCS Log replace by Id keyword, revision history cleaned up
	43	!
[1]	44	! Revision 1.25 2006/04/26 13:26:12 raasch
	45	! OpenMP optimization (+localsum, threadsum)
	46	!
	47	! Revision 1.1 1997/07/24 11:24:44 raasch
	48	! Initial revision
	49	!
	50	!
	51	! Description:
	52	! ------------
	53	! Compute the divergence of the provisional velocity field. Solve the Poisson
	54	! equation for the perturbation pressure. Compute the final velocities using
	55	! this perturbation pressure. Compute the remaining divergence.
	56	!------------------------------------------------------------------------------!
	57
	58	USE arrays_3d
	59	USE constants
	60	USE control_parameters
	61	USE cpulog
	62	USE grid_variables
	63	USE indices
	64	USE interfaces
	65	USE pegrid
	66	USE poisfft_mod
	67	USE poisfft_hybrid_mod
	68	USE statistics
	69
	70	IMPLICIT NONE
	71
	72	INTEGER :: i, j, k, sr
	73
[85]	74	REAL :: ddt_3d, localsum, threadsum
[1]	75
	76	REAL, DIMENSION(1:2) :: volume_flow_l, volume_flow_offset
[76]	77	REAL, DIMENSION(1:nzt) :: w_l, w_l_l
[1]	78
	79
	80	CALL cpu_log( log_point(8), 'pres', 'start' )
	81
[85]	82
	83	ddt_3d = 1.0 / dt_3d
	84
[1]	85	!
[667]	86	!-- Multigrid method expects 1 additional grid point for the arrays
	87	!-- d, tend and p
[1]	88	IF ( psolver == 'multigrid' ) THEN
[667]	89
[1]	90	DEALLOCATE( d )
[667]	91	ALLOCATE( d(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
	92
	93	IF ( ws_scheme_mom .OR. ws_scheme_sca ) THEN
	94
	95	DEALLOCATE( tend )
	96	ALLOCATE( tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
	97	DEALLOCATE( p )
	98	ALLOCATE( p(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
	99
	100	ENDIF
	101
[1]	102	ENDIF
	103
	104	!
[75]	105	!-- Conserve the volume flow at the outflow in case of non-cyclic lateral
	106	!-- boundary conditions
[106]	107	!-- WARNING: so far, this conservation does not work at the left/south
	108	!-- boundary if the topography at the inflow differs from that at the
	109	!-- outflow! For this case, volume_flow_area needs adjustment!
	110	!
	111	!-- Left/right
	112	IF ( conserve_volume_flow .AND. ( outflow_l .OR. outflow_r ) ) THEN
[75]	113
[106]	114	volume_flow(1) = 0.0
	115	volume_flow_l(1) = 0.0
	116
	117	IF ( outflow_l ) THEN
	118	i = 0
	119	ELSEIF ( outflow_r ) THEN
	120	i = nx+1
	121	ENDIF
	122
	123	DO j = nys, nyn
	124	!
	125	!-- Sum up the volume flow through the south/north boundary
	126	DO k = nzb_2d(j,i) + 1, nzt
[667]	127	volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k)
[106]	128	ENDDO
	129	ENDDO
	130
	131	#if defined( __parallel )
[622]	132	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[106]	133	CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 1, MPI_REAL, &
	134	MPI_SUM, comm1dy, ierr )
	135	#else
	136	volume_flow = volume_flow_l
	137	#endif
	138	volume_flow_offset(1) = ( volume_flow_initial(1) - volume_flow(1) ) &
	139	/ volume_flow_area(1)
	140
[667]	141	DO j = nysg, nyng
	142	DO k = nzb_2d(j,i) + 1, nzt
[106]	143	u(k,j,i) = u(k,j,i) + volume_flow_offset(1)
	144	ENDDO
	145	ENDDO
	146
	147	ENDIF
	148
	149	!
	150	!-- South/north
	151	IF ( conserve_volume_flow .AND. ( outflow_n .OR. outflow_s ) ) THEN
	152
[75]	153	volume_flow(2) = 0.0
	154	volume_flow_l(2) = 0.0
	155
[106]	156	IF ( outflow_s ) THEN
	157	j = 0
	158	ELSEIF ( outflow_n ) THEN
[75]	159	j = ny+1
[106]	160	ENDIF
	161
	162	DO i = nxl, nxr
[75]	163	!
[106]	164	!-- Sum up the volume flow through the south/north boundary
	165	DO k = nzb_2d(j,i) + 1, nzt
[667]	166	volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k)
[75]	167	ENDDO
[106]	168	ENDDO
	169
[75]	170	#if defined( __parallel )
[622]	171	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[75]	172	CALL MPI_ALLREDUCE( volume_flow_l(2), volume_flow(2), 1, MPI_REAL, &
	173	MPI_SUM, comm1dx, ierr )
	174	#else
	175	volume_flow = volume_flow_l
	176	#endif
	177	volume_flow_offset(2) = ( volume_flow_initial(2) - volume_flow(2) ) &
[106]	178	/ volume_flow_area(2)
[75]	179
[667]	180	DO i = nxlg, nxrg
[106]	181	DO k = nzb_v_inner(j,i) + 1, nzt
	182	v(k,j,i) = v(k,j,i) + volume_flow_offset(2)
[75]	183	ENDDO
[106]	184	ENDDO
[75]	185
	186	ENDIF
	187
[76]	188	!
	189	!-- Remove mean vertical velocity
	190	IF ( ibc_p_b == 1 .AND. ibc_p_t == 1 ) THEN
	191	IF ( simulated_time > 0.0 ) THEN ! otherwise nzb_w_inner is not yet known
	192	w_l = 0.0; w_l_l = 0.0
	193	DO i = nxl, nxr
	194	DO j = nys, nyn
	195	DO k = nzb_w_inner(j,i)+1, nzt
	196	w_l_l(k) = w_l_l(k) + w(k,j,i)
	197	ENDDO
	198	ENDDO
	199	ENDDO
	200	#if defined( __parallel )
[622]	201	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[76]	202	CALL MPI_ALLREDUCE( w_l_l(1), w_l(1), nzt, MPI_REAL, MPI_SUM, comm2d, &
	203	ierr )
	204	#else
	205	w_l = w_l_l
	206	#endif
	207	DO k = 1, nzt
	208	w_l(k) = w_l(k) / ngp_2dh_outer(k,0)
	209	ENDDO
[667]	210	DO i = nxlg, nxrg
	211	DO j = nysg, nyng
[76]	212	DO k = nzb_w_inner(j,i)+1, nzt
	213	w(k,j,i) = w(k,j,i) - w_l(k)
	214	ENDDO
	215	ENDDO
	216	ENDDO
	217	ENDIF
	218	ENDIF
[75]	219
	220	!
[1]	221	!-- Compute the divergence of the provisional velocity field.
	222	CALL cpu_log( log_point_s(1), 'divergence', 'start' )
	223
	224	IF ( psolver == 'multigrid' ) THEN
	225	!$OMP PARALLEL DO SCHEDULE( STATIC )
	226	DO i = nxl-1, nxr+1
	227	DO j = nys-1, nyn+1
	228	DO k = nzb, nzt+1
	229	d(k,j,i) = 0.0
	230	ENDDO
	231	ENDDO
	232	ENDDO
	233	ELSE
	234	!$OMP PARALLEL DO SCHEDULE( STATIC )
	235	DO i = nxl, nxra
	236	DO j = nys, nyna
	237	DO k = nzb+1, nzta
	238	d(k,j,i) = 0.0
	239	ENDDO
	240	ENDDO
	241	ENDDO
	242	ENDIF
	243
	244	localsum = 0.0
	245	threadsum = 0.0
	246
	247	#if defined( __ibm )
	248	!$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
	249	!$OMP DO SCHEDULE( STATIC )
	250	DO i = nxl, nxr
	251	DO j = nys, nyn
	252	DO k = nzb_s_inner(j,i)+1, nzt
[85]	253	d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
	254	( v(k,j+1,i) - v(k,j,i) ) * ddy + &
	255	( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ) * ddt_3d
[1]	256	ENDDO
	257	!
	258	!-- Additional pressure boundary condition at the bottom boundary for
	259	!-- inhomogeneous Prandtl layer heat fluxes and temperatures, respectively
	260	!-- dp/dz = -(dtau13/dx + dtau23/dy) + g*pt'/pt0.
	261	!-- This condition must not be applied at the start of a run, because then
	262	!-- flow_statistics has not yet been called and thus sums = 0.
	263	IF ( ibc_p_b == 2 .AND. sums(nzb+1,4) /= 0.0 ) THEN
	264	k = nzb_s_inner(j,i)
	265	d(k+1,j,i) = d(k+1,j,i) + ( &
	266	( usws(j,i+1) - usws(j,i) ) * ddx &
	267	+ ( vsws(j+1,i) - vsws(j,i) ) * ddy &
	268	- g * ( pt(k+1,j,i) - sums(k+1,4) ) / &
	269	sums(k+1,4) &
[85]	270	) * ddzw(k+1) * ddt_3d
[1]	271	ENDIF
	272
	273	!
	274	!-- Compute possible PE-sum of divergences for flow_statistics
	275	DO k = nzb_s_inner(j,i)+1, nzt
	276	threadsum = threadsum + ABS( d(k,j,i) )
	277	ENDDO
	278
	279	ENDDO
	280	ENDDO
	281
[85]	282	localsum = ( localsum + threadsum ) * dt_3d
[1]	283	!$OMP END PARALLEL
	284	#else
	285	IF ( ibc_p_b == 2 .AND. sums(nzb+1,4) /= 0.0 ) THEN
	286	!$OMP PARALLEL PRIVATE (i,j,k)
	287	!$OMP DO SCHEDULE( STATIC )
	288	DO i = nxl, nxr
	289	DO j = nys, nyn
	290	DO k = nzb_s_inner(j,i)+1, nzt
[667]	291	d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
	292	( v(k,j+1,i) - v(k,j,i) ) * ddy + &
	293	( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ) * ddt_3d
[1]	294	ENDDO
	295	ENDDO
	296	!
	297	!-- Additional pressure boundary condition at the bottom boundary for
	298	!-- inhomogeneous Prandtl layer heat fluxes and temperatures, respectively
	299	!-- dp/dz = -(dtau13/dx + dtau23/dy) + g*pt'/pt0.
	300	!-- This condition must not be applied at the start of a run, because then
	301	!-- flow_statistics has not yet been called and thus sums = 0.
	302	DO j = nys, nyn
	303	k = nzb_s_inner(j,i)
	304	d(k+1,j,i) = d(k+1,j,i) + ( &
	305	( usws(j,i+1) - usws(j,i) ) * ddx &
	306	+ ( vsws(j+1,i) - vsws(j,i) ) * ddy &
	307	- g * ( pt(k+1,j,i) - sums(k+1,4) ) / &
	308	sums(k+1,4) &
[85]	309	) * ddzw(k+1) * ddt_3d
[1]	310	ENDDO
	311	ENDDO
	312	!$OMP END PARALLEL
	313
	314	ELSE
	315
	316	!$OMP PARALLEL PRIVATE (i,j,k)
	317	!$OMP DO SCHEDULE( STATIC )
	318	DO i = nxl, nxr
	319	DO j = nys, nyn
	320	DO k = nzb_s_inner(j,i)+1, nzt
[85]	321	d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
[667]	322	( v(k,j+1,i) - v(k,j,i) ) * ddy + &
	323	( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ) * ddt_3d
[1]	324	ENDDO
	325	ENDDO
	326	ENDDO
	327	!$OMP END PARALLEL
	328
	329	ENDIF
	330
	331	!
	332	!-- Compute possible PE-sum of divergences for flow_statistics
	333	!$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
	334	!$OMP DO SCHEDULE( STATIC )
	335	DO i = nxl, nxr
	336	DO j = nys, nyn
	337	DO k = nzb+1, nzt
	338	threadsum = threadsum + ABS( d(k,j,i) )
	339	ENDDO
	340	ENDDO
	341	ENDDO
[85]	342	localsum = ( localsum + threadsum ) * dt_3d
[1]	343	!$OMP END PARALLEL
	344	#endif
	345
	346	!
	347	!-- For completeness, set the divergence sum of all statistic regions to those
	348	!-- of the total domain
	349	sums_divold_l(0:statistic_regions) = localsum
	350
	351	!
	352	!-- Determine absolute minimum/maximum (only for test cases, therefore as
	353	!-- comment line)
	354	! CALL global_min_max( nzb+1, nzt, nys, nyn, nxl, nxr, d, 'abs', divmax, &
[667]	355	! divmax_ijk )
[1]	356
	357	CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
	358
	359	!
	360	!-- Compute the pressure perturbation solving the Poisson equation
	361	IF ( psolver(1:7) == 'poisfft' ) THEN
	362
	363	!
	364	!-- Enlarge the size of tend, used as a working array for the transpositions
	365	IF ( nxra > nxr .OR. nyna > nyn .OR. nza > nz ) THEN
	366	DEALLOCATE( tend )
	367	ALLOCATE( tend(1:nza,nys:nyna,nxl:nxra) )
	368	ENDIF
	369
	370	!
	371	!-- Solve Poisson equation via FFT and solution of tridiagonal matrices
	372	IF ( psolver == 'poisfft' ) THEN
	373	!
	374	!-- Solver for 2d-decomposition
	375	CALL poisfft( d, tend )
	376	ELSEIF ( psolver == 'poisfft_hybrid' ) THEN
	377	!
	378	!-- Solver for 1d-decomposition (using MPI and OpenMP).
	379	!-- The old hybrid-solver is still included here, as long as there
	380	!-- are some optimization problems in poisfft
	381	CALL poisfft_hybrid( d )
	382	ENDIF
	383
	384	!
	385	!-- Resize tend to its normal size
	386	IF ( nxra > nxr .OR. nyna > nyn .OR. nza > nz ) THEN
	387	DEALLOCATE( tend )
[667]	388	ALLOCATE( tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]	389	ENDIF
	390
	391	!
	392	!-- Store computed perturbation pressure and set boundary condition in
	393	!-- z-direction
	394	!$OMP PARALLEL DO
	395	DO i = nxl, nxr
	396	DO j = nys, nyn
	397	DO k = nzb+1, nzt
	398	tend(k,j,i) = d(k,j,i)
	399	ENDDO
	400	ENDDO
	401	ENDDO
	402
	403	!
	404	!-- Bottom boundary:
	405	!-- This condition is only required for internal output. The pressure
	406	!-- gradient (dp(nzb+1)-dp(nzb))/dz is not used anywhere else.
	407	IF ( ibc_p_b == 1 ) THEN
	408	!
	409	!-- Neumann (dp/dz = 0)
	410	!$OMP PARALLEL DO
[667]	411	DO i = nxlg, nxrg
	412	DO j = nysg, nyng
[1]	413	tend(nzb_s_inner(j,i),j,i) = tend(nzb_s_inner(j,i)+1,j,i)
	414	ENDDO
	415	ENDDO
	416
	417	ELSEIF ( ibc_p_b == 2 ) THEN
	418	!
	419	!-- Neumann condition for inhomogeneous surfaces,
	420	!-- here currently still in the form of a zero gradient. Actually
	421	!-- dp/dz = -(dtau13/dx + dtau23/dy) + g*pt'/pt0 would have to be used for
	422	!-- the computation (cf. above: computation of divergences).
	423	!$OMP PARALLEL DO
[667]	424	DO i = nxlg, nxrg
	425	DO j = nysg, nyng
[1]	426	tend(nzb_s_inner(j,i),j,i) = tend(nzb_s_inner(j,i)+1,j,i)
	427	ENDDO
	428	ENDDO
	429
	430	ELSE
	431	!
	432	!-- Dirichlet
	433	!$OMP PARALLEL DO
[667]	434	DO i = nxlg, nxrg
	435	DO j = nysg, nyng
[1]	436	tend(nzb_s_inner(j,i),j,i) = 0.0
	437	ENDDO
	438	ENDDO
	439
	440	ENDIF
	441
	442	!
	443	!-- Top boundary
	444	IF ( ibc_p_t == 1 ) THEN
	445	!
	446	!-- Neumann
	447	!$OMP PARALLEL DO
[667]	448	DO i = nxlg, nxrg
	449	DO j = nysg, nyng
[1]	450	tend(nzt+1,j,i) = tend(nzt,j,i)
	451	ENDDO
	452	ENDDO
	453
	454	ELSE
	455	!
	456	!-- Dirichlet
	457	!$OMP PARALLEL DO
[667]	458	DO i = nxlg, nxrg
	459	DO j = nysg, nyng
[1]	460	tend(nzt+1,j,i) = 0.0
	461	ENDDO
	462	ENDDO
	463
	464	ENDIF
	465
	466	!
	467	!-- Exchange boundaries for p
[667]	468	CALL exchange_horiz( tend, nbgp )
[1]	469
	470	ELSEIF ( psolver == 'sor' ) THEN
	471
	472	!
	473	!-- Solve Poisson equation for perturbation pressure using SOR-Red/Black
	474	!-- scheme
[667]	475	CALL sor( d, ddzu_pres, ddzw, p )
[1]	476	tend = p
	477
	478	ELSEIF ( psolver == 'multigrid' ) THEN
	479
	480	!
	481	!-- Solve Poisson equation for perturbation pressure using Multigrid scheme,
[667]	482	!-- array tend is used to store the residuals, logical exchange_mg is used
	483	!-- to discern data exchange in multigrid ( 1 ghostpoint ) and normal grid
	484	!-- ( nbgp ghost points ).
	485	exchange_mg = .TRUE.
[1]	486	CALL poismg( tend )
[667]	487	exchange_mg = .FALSE.
[1]	488	!
	489	!-- Restore perturbation pressure on tend because this array is used
	490	!-- further below to correct the velocity fields
[667]	491
[1]	492	tend = p
[667]	493	IF( ws_scheme_mom .OR. ws_scheme_sca ) THEN
	494	!
	495	!-- Allocate p to its normal size and restore pressure.
	496	DEALLOCATE( p )
	497	ALLOCATE( p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	498	DO i = nxl, nxr
	499	DO j = nys, nyn
	500	DO k = nzb_s_inner(j,i), nzt
	501	p(k,j,i) = tend(k,j,i)
	502	ENDDO
	503	ENDDO
	504	ENDDO
	505	ENDIF
[1]	506
	507	ENDIF
	508
	509	!
	510	!-- Store perturbation pressure on array p, used in the momentum equations
	511	IF ( psolver(1:7) == 'poisfft' ) THEN
	512	!
	513	!-- Here, only the values from the left and right boundaries are copied
	514	!-- The remaining values are copied in the following loop due to speed
	515	!-- optimization
	516	!$OMP PARALLEL DO
[667]	517	DO j = nysg, nyng
[1]	518	DO k = nzb, nzt+1
[667]	519	p(k,j,nxlg:nxl-1) = tend(k,j,nxlg:nxl-1)
	520	p(k,j,nxr+1:nxrg) = tend(k,j,nxr+1:nxrg)
[1]	521	ENDDO
	522	ENDDO
	523	ENDIF
	524
	525	!
	526	!-- Correction of the provisional velocities with the current perturbation
	527	!-- pressure just computed
[75]	528	IF ( conserve_volume_flow .AND. &
	529	( bc_lr == 'cyclic' .OR. bc_ns == 'cyclic' ) ) THEN
[1]	530	volume_flow_l(1) = 0.0
	531	volume_flow_l(2) = 0.0
	532	ENDIF
	533	!$OMP PARALLEL PRIVATE (i,j,k)
	534	!$OMP DO
	535	DO i = nxl, nxr
	536	IF ( psolver(1:7) == 'poisfft' ) THEN
[667]	537	DO j = nysg, nyng
[1]	538	DO k = nzb, nzt+1
	539	p(k,j,i) = tend(k,j,i)
	540	ENDDO
	541	ENDDO
	542	ENDIF
	543	DO j = nys, nyn
	544	DO k = nzb_w_inner(j,i)+1, nzt
	545	w(k,j,i) = w(k,j,i) - dt_3d * &
	546	( tend(k+1,j,i) - tend(k,j,i) ) * ddzu(k+1)
	547	ENDDO
	548	DO k = nzb_u_inner(j,i)+1, nzt
	549	u(k,j,i) = u(k,j,i) - dt_3d * ( tend(k,j,i) - tend(k,j,i-1) ) * ddx
	550	ENDDO
	551	DO k = nzb_v_inner(j,i)+1, nzt
	552	v(k,j,i) = v(k,j,i) - dt_3d * ( tend(k,j,i) - tend(k,j-1,i) ) * ddy
	553	ENDDO
	554
	555	!
	556	!-- Sum up the volume flow through the right and north boundary
[75]	557	IF ( conserve_volume_flow .AND. bc_lr == 'cyclic' .AND. &
[667]	558	bc_ns == 'cyclic' .AND. i == nx ) THEN
[1]	559	!$OMP CRITICAL
	560	DO k = nzb_2d(j,i) + 1, nzt
[667]	561	volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k)
[1]	562	ENDDO
	563	!$OMP END CRITICAL
	564	ENDIF
[75]	565	IF ( conserve_volume_flow .AND. bc_ns == 'cyclic' .AND. &
[667]	566	bc_lr == 'cyclic' .AND. j == ny ) THEN
[1]	567	!$OMP CRITICAL
	568	DO k = nzb_2d(j,i) + 1, nzt
[667]	569	volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k)
[1]	570	ENDDO
	571	!$OMP END CRITICAL
	572	ENDIF
	573
	574	ENDDO
	575	ENDDO
	576	!$OMP END PARALLEL
	577
	578	!
[667]	579	!-- Resize tend to its normal size in case of multigrid and ws-scheme.
	580	IF ( psolver == 'multigrid' .AND. ( ws_scheme_mom &
	581	.OR. ws_scheme_sca ) ) THEN
	582	DEALLOCATE( tend )
	583	ALLOCATE( tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	584	ENDIF
	585
	586	!
[1]	587	!-- Conserve the volume flow
[75]	588	IF ( conserve_volume_flow .AND. &
[667]	589	( bc_lr == 'cyclic' .AND. bc_ns == 'cyclic' ) ) THEN
[1]	590
	591	#if defined( __parallel )
[622]	592	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[1]	593	CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 2, MPI_REAL, &
	594	MPI_SUM, comm2d, ierr )
	595	#else
	596	volume_flow = volume_flow_l
	597	#endif
	598
	599	volume_flow_offset = ( volume_flow_initial - volume_flow ) / &
	600	volume_flow_area
	601
	602	!$OMP PARALLEL PRIVATE (i,j,k)
	603	!$OMP DO
	604	DO i = nxl, nxr
	605	DO j = nys, nyn
[667]	606	DO k = nzb_u_inner(j,i) + 1, nzt
	607	u(k,j,i) = u(k,j,i) + volume_flow_offset(1)
	608	v(k,j,i) = v(k,j,i) + volume_flow_offset(2)
	609	ENDDO
[1]	610	ENDDO
	611	ENDDO
[667]	612
[1]	613	!$OMP END PARALLEL
	614
	615	ENDIF
	616
	617	!
	618	!-- Exchange of boundaries for the velocities
[667]	619	CALL exchange_horiz( u, nbgp )
	620	CALL exchange_horiz( v, nbgp )
	621	CALL exchange_horiz( w, nbgp )
[1]	622
	623	!
	624	!-- Compute the divergence of the corrected velocity field,
	625	!-- a possible PE-sum is computed in flow_statistics
	626	CALL cpu_log( log_point_s(1), 'divergence', 'start' )
	627	sums_divnew_l = 0.0
	628
	629	!
	630	!-- d must be reset to zero because it can contain nonzero values below the
	631	!-- topography
	632	IF ( topography /= 'flat' ) d = 0.0
	633
	634	localsum = 0.0
	635	threadsum = 0.0
	636
	637	!$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
	638	!$OMP DO SCHEDULE( STATIC )
	639	#if defined( __ibm )
	640	DO i = nxl, nxr
	641	DO j = nys, nyn
	642	DO k = nzb_s_inner(j,i)+1, nzt
	643	d(k,j,i) = ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
	644	( v(k,j+1,i) - v(k,j,i) ) * ddy + &
	645	( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	646	ENDDO
	647	DO k = nzb+1, nzt
	648	threadsum = threadsum + ABS( d(k,j,i) )
	649	ENDDO
	650	ENDDO
	651	ENDDO
	652	#else
	653	DO i = nxl, nxr
	654	DO j = nys, nyn
	655	DO k = nzb_s_inner(j,i)+1, nzt
	656	d(k,j,i) = ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
	657	( v(k,j+1,i) - v(k,j,i) ) * ddy + &
	658	( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	659	threadsum = threadsum + ABS( d(k,j,i) )
	660	ENDDO
	661	ENDDO
	662	ENDDO
	663	#endif
[667]	664
[1]	665	localsum = localsum + threadsum
	666	!$OMP END PARALLEL
	667
	668	!
	669	!-- For completeness, set the divergence sum of all statistic regions to those
	670	!-- of the total domain
	671	sums_divnew_l(0:statistic_regions) = localsum
	672
	673	CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
	674
	675	CALL cpu_log( log_point(8), 'pres', 'stop' )
[667]	676
[1]	677
	678
	679	END SUBROUTINE pres

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |