Home

Context Navigation

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

Last change on this file since 692 was 682, checked in by gryschka, 14 years ago
removed print-statement for debugging
Property svn:keywords set to `Id`
File size: 23.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |