Home

Context Navigation

pres.f90 @ 622

Last change on this file since 622 was 622, checked in by raasch, 13 years ago

New:
---

Optional barriers included in order to speed up collective operations
MPI_ALLTOALL and MPI_ALLREDUCE. This feature is controlled with new initial
parameter collective_wait. Default is .FALSE, but .TRUE. on SGI-type
systems. (advec_particles, advec_s_bc, buoyancy, check_for_restart,
cpu_statistics, data_output_2d, data_output_ptseries, flow_statistics,
global_min_max, inflow_turbulence, init_3d_model, init_particles, init_pegrid,
init_slope, parin, pres, poismg, set_particle_attributes, timestep,
read_var_list, user_statistics, write_compressed, write_var_list)

Adjustments for Kyushu Univ. (lcrte, ibmku). Concerning hybrid
(MPI/openMP) runs, the number of openMP threads per MPI tasks can now
be given as an argument to mrun-option -O. (mbuild, mrun, subjob)

Changed:

Initialization of the module command changed for SGI-ICE/lcsgi (mbuild, subjob)

Errors:

Property svn:keywords set to Id

File size: 19.0 KB

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

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

Context Navigation

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

Download in other formats: