Home

Context Navigation

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

Last change on this file since 698 was 695, checked in by raasch, 14 years ago
comment line changed
Property svn:keywords set to `Id`
File size: 23.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |