Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |