Home

Context Navigation

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

Last change on this file since 3401 was 3347, checked in by suehring, 6 years ago
Offline nesting revised and separated from large_scale_forcing_mod; Time-dependent synthetic turbulence generator; bugfixes in USM/LSM radiation model and init_pegrid
Property svn:keywords set to `Id`
File size: 32.2 KB

Rev	Line
[1682]	1	!> @file pres.f90
[2000]	2	!------------------------------------------------------------------------------!
[2696]	3	! This file is part of the PALM model system.
[1036]	4	!
[2000]	5	! PALM is free software: you can redistribute it and/or modify it under the
	6	! terms of the GNU General Public License as published by the Free Software
	7	! Foundation, either version 3 of the License, or (at your option) any later
	8	! version.
[1036]	9	!
	10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
	11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
	12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
	13	!
	14	! You should have received a copy of the GNU General Public License along with
	15	! PALM. If not, see <http://www.gnu.org/licenses/>.
	16	!
[2718]	17	! Copyright 1997-2018 Leibniz Universitaet Hannover
[2000]	18	!------------------------------------------------------------------------------!
[1036]	19	!
[484]	20	! Current revisions:
[1212]	21	! ------------------
[1932]	22	!
[3183]	23	!
[1919]	24	! Former revisions:
	25	! -----------------
[3016]	26	! $Id: pres.f90 3347 2018-10-15 14:21:08Z eckhard $
[3347]	27	! Bugfixes in offline nesting.
	28	! Add comment.
	29	!
	30	! 3341 2018-10-15 10:31:27Z suehring
[3183]	31	! Rename variables for boundary flags and nesting
	32	!
	33	! 3182 2018-07-27 13:36:03Z suehring
[3016]	34	! Dollar sign added before Id
	35	!
	36	! 2696 2017-12-14 17:12:51Z kanani
[3013]	37	! To avoid jumps while plotting w-profiles w level nzt+1 is set to w level nzt
	38	! after velocity modifications through the pressure solver were carried out
	39	!
	40	! 2696 2017-12-14 17:12:51Z kanani
[2716]	41	! Corrected "Former revisions" section
[2696]	42	!
[2716]	43	! 2696 2017-12-14 17:12:51Z kanani
	44	! Change in file header (GPL part)
	45	! poismg_noopt modularized (MS)
	46	!
[2696]	47	! 2298 2017-06-29 09:28:18Z raasch
[2298]	48	! comment changed + some formatting
	49	!
	50	! 2233 2017-05-30 18:08:54Z suehring
[1919]	51	!
[2233]	52	! 2232 2017-05-30 17:47:52Z suehring
	53	! Adjustments to new topography and surface concept
	54	!
[2119]	55	! 2118 2017-01-17 16:38:49Z raasch
	56	! OpenACC directives and related code removed
	57	!
[2074]	58	! 2073 2016-11-30 14:34:05Z raasch
	59	! openmp bugfix for calculation of new divergence
	60	!
[2038]	61	! 2037 2016-10-26 11:15:40Z knoop
	62	! Anelastic approximation implemented
	63	!
[2001]	64	! 2000 2016-08-20 18:09:15Z knoop
	65	! Forced header and separation lines into 80 columns
	66	!
[1933]	67	! 1932 2016-06-10 12:09:21Z suehring
	68	! Initial version of purely vertical nesting introduced.
	69	!
[1932]	70	! 1931 2016-06-10 12:06:59Z suehring
	71	! Rename multigrid into multigrid_noopt and multigrid_fast into multigrid
	72	!
[1930]	73	! 1929 2016-06-09 16:25:25Z suehring
	74	! Bugfix: weight_substep for initial call, replace by local variable
	75	!
[1919]	76	! 1918 2016-05-27 14:35:57Z raasch
[1918]	77	! sum of divergence is also calculated when pres is called before the initial
	78	! first time step,
	79	! stearing is modified by using intermediate_timestep_count = 0 in order to
	80	! determine if pres is called before the first initial timestep,
	81	! bugfix: in case of Neumann conditions for pressure both at bottom and top,
	82	! mean vertical velocity is also removed for the first time step
	83	! bugfix for calculating divergences
[1321]	84	!
[1909]	85	! 1908 2016-05-25 17:22:32Z suehring
	86	! New divergence for output into RUN_CONTROL file is calculated only at last
	87	! Runge-Kutta step
	88	!
[1846]	89	! 1845 2016-04-08 08:29:13Z raasch
	90	! nzb_2d replace by nzb_u\|v_inner
	91	!
[1800]	92	! 1799 2016-04-05 08:35:55Z gronemeier
	93	! Bugfix: excluded third dimension from horizontal volume flow calculation
	94	!
[1763]	95	! 1762 2016-02-25 12:31:13Z hellstea
	96	! Introduction of nested domain feature
	97	!
[1683]	98	! 1682 2015-10-07 23:56:08Z knoop
	99	! Code annotations made doxygen readable
	100	!
[1576]	101	! 1575 2015-03-27 09:56:27Z raasch
	102	! poismg_fast + respective module added, adjustments for psolver-queries
	103	!
[1343]	104	! 1342 2014-03-26 17:04:47Z kanani
	105	! REAL constants defined as wp-kind
	106	!
[1321]	107	! 1320 2014-03-20 08:40:49Z raasch
[1320]	108	! ONLY-attribute added to USE-statements,
	109	! kind-parameters added to all INTEGER and REAL declaration statements,
	110	! kinds are defined in new module kinds,
	111	! old module precision_kind is removed,
	112	! revision history before 2012 removed,
	113	! comment fields (!:) to be used for variable explanations added to
	114	! all variable declaration statements
[708]	115	!
[1319]	116	! 1318 2014-03-17 13:35:16Z raasch
	117	! module interfaces removed
	118	!
[1307]	119	! 1306 2014-03-13 14:30:59Z raasch
	120	! second argument removed from routine poisfft
	121	!
[1258]	122	! 1257 2013-11-08 15:18:40Z raasch
	123	! openacc loop and loop vector clauses removed, independent clauses added,
	124	! end parallel replaced by end parallel loop
	125	!
[1222]	126	! 1221 2013-09-10 08:59:13Z raasch
	127	! openACC porting of reduction operations, loops for calculating d are
	128	! using the rflags_s_inner multiply flag instead of the nzb_s_inner loop index
	129	!
[1213]	130	! 1212 2013-08-15 08:46:27Z raasch
	131	! call of poisfft_hybrid removed
	132	!
[1118]	133	! 1117 2013-03-27 11:15:36Z suehring
	134	! Bugfix in OpenMP parallelization.
	135	!
[1114]	136	! 1113 2013-03-10 02:48:14Z raasch
	137	! GPU-porting of several loops, some loops rearranged
	138	!
[1112]	139	! 1111 2013-03-08 23:54:10Z
	140	! openACC statements added,
	141	! ibc_p_b = 2 removed
	142	!
[1093]	143	! 1092 2013-02-02 11:24:22Z raasch
	144	! unused variables removed
	145	!
[1037]	146	! 1036 2012-10-22 13:43:42Z raasch
	147	! code put under GPL (PALM 3.9)
	148	!
[1004]	149	! 1003 2012-09-14 14:35:53Z raasch
	150	! adjustment of array tend for cases with unequal subdomain sizes removed
	151	!
[1]	152	! Revision 1.1 1997/07/24 11:24:44 raasch
	153	! Initial revision
	154	!
	155	!
	156	! Description:
	157	! ------------
[1682]	158	!> Compute the divergence of the provisional velocity field. Solve the Poisson
	159	!> equation for the perturbation pressure. Compute the final velocities using
	160	!> this perturbation pressure. Compute the remaining divergence.
[1]	161	!------------------------------------------------------------------------------!
[1682]	162	SUBROUTINE pres
	163
[1]	164
[1320]	165	USE arrays_3d, &
[2037]	166	ONLY: d, ddzu, ddzu_pres, ddzw, dzw, p, p_loc, rho_air, rho_air_zw, &
	167	tend, u, v, w
[1320]	168
	169	USE control_parameters, &
[3182]	170	ONLY: bc_lr_cyc, bc_ns_cyc, bc_radiation_l, bc_radiation_n, &
	171	bc_radiation_r, bc_radiation_s, child_domain, &
	172	conserve_volume_flow, coupling_mode, &
[3012]	173	dt_3d, gathered_size, ibc_p_b, ibc_p_t, &
	174	intermediate_timestep_count, intermediate_timestep_count_max, &
[3347]	175	mg_switch_to_pe0_level, nesting_offline, &
	176	psolver, subdomain_size, &
[3182]	177	topography, volume_flow, volume_flow_area, volume_flow_initial
[1320]	178
	179	USE cpulog, &
	180	ONLY: cpu_log, log_point, log_point_s
	181
	182	USE grid_variables, &
	183	ONLY: ddx, ddy
	184
	185	USE indices, &
	186	ONLY: nbgp, ngp_2dh_outer, nx, nxl, nxlg, nxl_mg, nxr, nxrg, nxr_mg, &
[2232]	187	ny, nys, nysg, nys_mg, nyn, nyng, nyn_mg, nzb, nzt, nzt_mg, &
	188	wall_flags_0
[1320]	189
	190	USE kinds
	191
[1]	192	USE pegrid
[1933]	193
	194	USE pmc_interface, &
	195	ONLY: nesting_mode
[1]	196
[1320]	197	USE poisfft_mod, &
	198	ONLY: poisfft
	199
[1575]	200	USE poismg_mod
	201
[2696]	202	USE poismg_noopt_mod
	203
[1320]	204	USE statistics, &
	205	ONLY: statistic_regions, sums_divnew_l, sums_divold_l, weight_pres, &
	206	weight_substep
	207
[2232]	208	USE surface_mod, &
	209	ONLY : bc_h
	210
[1]	211	IMPLICIT NONE
	212
[1682]	213	INTEGER(iwp) :: i !<
	214	INTEGER(iwp) :: j !<
	215	INTEGER(iwp) :: k !<
[2232]	216	INTEGER(iwp) :: m !<
[1]	217
[1682]	218	REAL(wp) :: ddt_3d !<
[1918]	219	REAL(wp) :: d_weight_pres !<
[1682]	220	REAL(wp) :: localsum !<
	221	REAL(wp) :: threadsum !<
[1918]	222	REAL(wp) :: weight_pres_l !<
[1929]	223	REAL(wp) :: weight_substep_l !<
[1]	224
[1762]	225	REAL(wp), DIMENSION(1:3) :: volume_flow_l !<
	226	REAL(wp), DIMENSION(1:3) :: volume_flow_offset !<
[1682]	227	REAL(wp), DIMENSION(1:nzt) :: w_l !<
	228	REAL(wp), DIMENSION(1:nzt) :: w_l_l !<
[1]	229
[3182]	230	LOGICAL :: child_domain_nvn !<
[1]	231
[1933]	232
[1]	233	CALL cpu_log( log_point(8), 'pres', 'start' )
	234
[85]	235
[1918]	236	!
	237	!-- Calculate quantities to be used locally
[1342]	238	ddt_3d = 1.0_wp / dt_3d
[1918]	239	IF ( intermediate_timestep_count == 0 ) THEN
	240	!
	241	!-- If pres is called before initial time step
[1929]	242	weight_pres_l = 1.0_wp
	243	d_weight_pres = 1.0_wp
	244	weight_substep_l = 1.0_wp
[1918]	245	ELSE
[1929]	246	weight_pres_l = weight_pres(intermediate_timestep_count)
	247	d_weight_pres = 1.0_wp / weight_pres(intermediate_timestep_count)
	248	weight_substep_l = weight_substep(intermediate_timestep_count)
[1918]	249	ENDIF
[85]	250
[1]	251	!
[707]	252	!-- Multigrid method expects array d to have one ghost layer.
	253	!--
[1575]	254	IF ( psolver(1:9) == 'multigrid' ) THEN
[667]	255
[1]	256	DEALLOCATE( d )
[667]	257	ALLOCATE( d(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
[707]	258
	259	!
	260	!-- Since p is later used to hold the weighted average of the substeps, it
	261	!-- cannot be used in the iterative solver. Therefore, its initial value is
	262	!-- stored on p_loc, which is then iteratively advanced in every substep.
[1762]	263	IF ( intermediate_timestep_count <= 1 ) THEN
[707]	264	DO i = nxl-1, nxr+1
	265	DO j = nys-1, nyn+1
	266	DO k = nzb, nzt+1
	267	p_loc(k,j,i) = p(k,j,i)
	268	ENDDO
	269	ENDDO
	270	ENDDO
[667]	271	ENDIF
	272
[1918]	273	ELSEIF ( psolver == 'sor' .AND. intermediate_timestep_count <= 1 ) THEN
[707]	274
	275	!
	276	!-- Since p is later used to hold the weighted average of the substeps, it
	277	!-- cannot be used in the iterative solver. Therefore, its initial value is
	278	!-- stored on p_loc, which is then iteratively advanced in every substep.
	279	p_loc = p
	280
[1]	281	ENDIF
	282
	283	!
[75]	284	!-- Conserve the volume flow at the outflow in case of non-cyclic lateral
	285	!-- boundary conditions
[106]	286	!-- WARNING: so far, this conservation does not work at the left/south
	287	!-- boundary if the topography at the inflow differs from that at the
	288	!-- outflow! For this case, volume_flow_area needs adjustment!
	289	!
	290	!-- Left/right
[3182]	291	IF ( conserve_volume_flow .AND. ( bc_radiation_l .OR. &
	292	bc_radiation_r ) ) THEN
[680]	293
[1342]	294	volume_flow(1) = 0.0_wp
	295	volume_flow_l(1) = 0.0_wp
[106]	296
[3182]	297	IF ( bc_radiation_l ) THEN
[106]	298	i = 0
[3182]	299	ELSEIF ( bc_radiation_r ) THEN
[106]	300	i = nx+1
	301	ENDIF
	302
	303	DO j = nys, nyn
	304	!
	305	!-- Sum up the volume flow through the south/north boundary
[2232]	306	DO k = nzb+1, nzt
	307	volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k) &
	308	* MERGE( 1.0_wp, 0.0_wp, &
	309	BTEST( wall_flags_0(k,j,i), 1 ) &
	310	)
[106]	311	ENDDO
	312	ENDDO
	313
	314	#if defined( __parallel )
[680]	315	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
[106]	316	CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 1, MPI_REAL, &
	317	MPI_SUM, comm1dy, ierr )
	318	#else
	319	volume_flow = volume_flow_l
	320	#endif
[709]	321	volume_flow_offset(1) = ( volume_flow_initial(1) - volume_flow(1) ) &
[106]	322	/ volume_flow_area(1)
	323
[667]	324	DO j = nysg, nyng
[2232]	325	DO k = nzb+1, nzt
	326	u(k,j,i) = u(k,j,i) + volume_flow_offset(1) &
	327	* MERGE( 1.0_wp, 0.0_wp, &
	328	BTEST( wall_flags_0(k,j,i), 1 ) &
	329	)
[106]	330	ENDDO
	331	ENDDO
	332
	333	ENDIF
	334
	335	!
	336	!-- South/north
[3182]	337	IF ( conserve_volume_flow .AND. ( bc_radiation_n .OR. bc_radiation_s ) ) THEN
[106]	338
[1342]	339	volume_flow(2) = 0.0_wp
	340	volume_flow_l(2) = 0.0_wp
[75]	341
[3182]	342	IF ( bc_radiation_s ) THEN
[106]	343	j = 0
[3182]	344	ELSEIF ( bc_radiation_n ) THEN
[75]	345	j = ny+1
[106]	346	ENDIF
	347
	348	DO i = nxl, nxr
[75]	349	!
[106]	350	!-- Sum up the volume flow through the south/north boundary
[2232]	351	DO k = nzb+1, nzt
	352	volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k) &
	353	* MERGE( 1.0_wp, 0.0_wp, &
	354	BTEST( wall_flags_0(k,j,i), 2 ) &
	355	)
[75]	356	ENDDO
[106]	357	ENDDO
	358
[75]	359	#if defined( __parallel )
[680]	360	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
[75]	361	CALL MPI_ALLREDUCE( volume_flow_l(2), volume_flow(2), 1, MPI_REAL, &
	362	MPI_SUM, comm1dx, ierr )
	363	#else
	364	volume_flow = volume_flow_l
	365	#endif
	366	volume_flow_offset(2) = ( volume_flow_initial(2) - volume_flow(2) ) &
[106]	367	/ volume_flow_area(2)
[75]	368
[667]	369	DO i = nxlg, nxrg
[2232]	370	DO k = nzb+1, nzt
	371	v(k,j,i) = v(k,j,i) + volume_flow_offset(2) &
	372	* MERGE( 1.0_wp, 0.0_wp, &
	373	BTEST( wall_flags_0(k,j,i), 2 ) &
	374	)
[75]	375	ENDDO
[106]	376	ENDDO
[75]	377
	378	ENDIF
	379
[76]	380	!
[1762]	381	!-- Remove mean vertical velocity in case that Neumann conditions are
[1933]	382	!-- used both at bottom and top boundary, and if not a nested domain in a
	383	!-- normal nesting run. In case of vertical nesting, this must be done.
[3182]	384	!-- Therefore an auxiliary logical variable child_domain_nvn is used here, and
[1933]	385	!-- nvn stands for non-vertical nesting.
[1918]	386	!-- This cannot be done before the first initial time step because ngp_2dh_outer
	387	!-- is not yet known then.
[3182]	388	child_domain_nvn = child_domain
	389	IF ( child_domain .AND. nesting_mode == 'vertical' ) THEN
	390	child_domain_nvn = .FALSE.
[1933]	391	ENDIF
	392
[3347]	393	IF ( ibc_p_b == 1 .AND. ibc_p_t == 1 .AND. .NOT. nesting_offline .AND. &
[3182]	394	.NOT. child_domain_nvn .AND. intermediate_timestep_count /= 0 ) &
[1918]	395	THEN
	396	w_l = 0.0_wp; w_l_l = 0.0_wp
	397	DO i = nxl, nxr
	398	DO j = nys, nyn
[2232]	399	DO k = nzb+1, nzt
	400	w_l_l(k) = w_l_l(k) + w(k,j,i) &
	401	* MERGE( 1.0_wp, 0.0_wp, &
	402	BTEST( wall_flags_0(k,j,i), 3 ) &
	403	)
[76]	404	ENDDO
	405	ENDDO
[1918]	406	ENDDO
[76]	407	#if defined( __parallel )
[1918]	408	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
	409	CALL MPI_ALLREDUCE( w_l_l(1), w_l(1), nzt, MPI_REAL, MPI_SUM, &
	410	comm2d, ierr )
[76]	411	#else
[1918]	412	w_l = w_l_l
[76]	413	#endif
[1918]	414	DO k = 1, nzt
	415	w_l(k) = w_l(k) / ngp_2dh_outer(k,0)
	416	ENDDO
	417	DO i = nxlg, nxrg
	418	DO j = nysg, nyng
[2232]	419	DO k = nzb+1, nzt
	420	w(k,j,i) = w(k,j,i) - w_l(k) &
	421	* MERGE( 1.0_wp, 0.0_wp, &
	422	BTEST( wall_flags_0(k,j,i), 3 ) &
	423	)
[76]	424	ENDDO
	425	ENDDO
[1918]	426	ENDDO
[76]	427	ENDIF
[75]	428
	429	!
[1]	430	!-- Compute the divergence of the provisional velocity field.
	431	CALL cpu_log( log_point_s(1), 'divergence', 'start' )
	432
[1575]	433	IF ( psolver(1:9) == 'multigrid' ) THEN
[2232]	434	!$OMP PARALLEL DO SCHEDULE( STATIC ) PRIVATE (i,j,k)
[1]	435	DO i = nxl-1, nxr+1
	436	DO j = nys-1, nyn+1
	437	DO k = nzb, nzt+1
[1342]	438	d(k,j,i) = 0.0_wp
[1]	439	ENDDO
	440	ENDDO
	441	ENDDO
	442	ELSE
[2232]	443	!$OMP PARALLEL DO SCHEDULE( STATIC ) PRIVATE (i,j,k)
[1003]	444	DO i = nxl, nxr
	445	DO j = nys, nyn
	446	DO k = nzb+1, nzt
[1342]	447	d(k,j,i) = 0.0_wp
[1]	448	ENDDO
	449	ENDDO
	450	ENDDO
	451	ENDIF
	452
[1342]	453	localsum = 0.0_wp
	454	threadsum = 0.0_wp
[1]	455
	456	#if defined( __ibm )
	457	!$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
	458	!$OMP DO SCHEDULE( STATIC )
	459	DO i = nxl, nxr
	460	DO j = nys, nyn
[2232]	461	DO k = nzb+1, nzt
[2037]	462	d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx + &
	463	( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy + &
	464	( w(k,j,i) * rho_air_zw(k) - &
	465	w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k) &
[2232]	466	) * ddt_3d * d_weight_pres &
	467	* MERGE( 1.0_wp, 0.0_wp, &
	468	BTEST( wall_flags_0(k,j,i), 0 ) &
	469	)
[1]	470	ENDDO
	471	!
	472	!-- Compute possible PE-sum of divergences for flow_statistics
[2232]	473	DO k = nzb+1, nzt
	474	threadsum = threadsum + ABS( d(k,j,i) ) &
	475	* MERGE( 1.0_wp, 0.0_wp, &
	476	BTEST( wall_flags_0(k,j,i), 0 ) &
	477	)
[1]	478	ENDDO
	479
	480	ENDDO
	481	ENDDO
	482
[1918]	483	IF ( intermediate_timestep_count == intermediate_timestep_count_max .OR. &
	484	intermediate_timestep_count == 0 ) THEN
	485	localsum = localsum + threadsum * dt_3d * weight_pres_l
[1908]	486	ENDIF
[1]	487	!$OMP END PARALLEL
	488	#else
	489
[1111]	490	!$OMP PARALLEL PRIVATE (i,j,k)
	491	!$OMP DO SCHEDULE( STATIC )
	492	DO i = nxl, nxr
	493	DO j = nys, nyn
	494	DO k = 1, nzt
[2037]	495	d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx + &
	496	( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy + &
	497	( w(k,j,i) * rho_air_zw(k) - &
	498	w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k) &
[2232]	499	) * ddt_3d * d_weight_pres &
	500	* MERGE( 1.0_wp, 0.0_wp, &
	501	BTEST( wall_flags_0(k,j,i), 0 ) &
	502	)
[1]	503	ENDDO
	504	ENDDO
[1111]	505	ENDDO
	506	!$OMP END PARALLEL
[1]	507
	508	!
[1908]	509	!-- Compute possible PE-sum of divergences for flow_statistics. Carry out
	510	!-- computation only at last Runge-Kutta substep.
[1918]	511	IF ( intermediate_timestep_count == intermediate_timestep_count_max .OR. &
	512	intermediate_timestep_count == 0 ) THEN
[1908]	513	!$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
	514	!$OMP DO SCHEDULE( STATIC )
	515	DO i = nxl, nxr
	516	DO j = nys, nyn
	517	DO k = nzb+1, nzt
	518	threadsum = threadsum + ABS( d(k,j,i) )
	519	ENDDO
[1]	520	ENDDO
	521	ENDDO
[1918]	522	localsum = localsum + threadsum * dt_3d * weight_pres_l
[1908]	523	!$OMP END PARALLEL
	524	ENDIF
[1]	525	#endif
	526
	527	!
	528	!-- For completeness, set the divergence sum of all statistic regions to those
	529	!-- of the total domain
[1918]	530	IF ( intermediate_timestep_count == intermediate_timestep_count_max .OR. &
	531	intermediate_timestep_count == 0 ) THEN
[1908]	532	sums_divold_l(0:statistic_regions) = localsum
[1918]	533	ENDIF
[1]	534
	535	CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
	536
	537	!
	538	!-- Compute the pressure perturbation solving the Poisson equation
[1575]	539	IF ( psolver == 'poisfft' ) THEN
[1]	540
	541	!
	542	!-- Solve Poisson equation via FFT and solution of tridiagonal matrices
[1575]	543	CALL poisfft( d )
[1212]	544
[1]	545	!
	546	!-- Store computed perturbation pressure and set boundary condition in
	547	!-- z-direction
[2232]	548	!$OMP PARALLEL DO PRIVATE (i,j,k)
[1]	549	DO i = nxl, nxr
	550	DO j = nys, nyn
	551	DO k = nzb+1, nzt
	552	tend(k,j,i) = d(k,j,i)
	553	ENDDO
	554	ENDDO
	555	ENDDO
	556
	557	!
	558	!-- Bottom boundary:
	559	!-- This condition is only required for internal output. The pressure
	560	!-- gradient (dp(nzb+1)-dp(nzb))/dz is not used anywhere else.
	561	IF ( ibc_p_b == 1 ) THEN
	562	!
[2232]	563	!-- Neumann (dp/dz = 0). Using surfae data type, first for non-natural
	564	!-- surfaces, then for natural and urban surfaces
	565	!-- Upward facing
	566	!$OMP PARALLEL DO PRIVATE( i, j, k )
	567	DO m = 1, bc_h(0)%ns
[2298]	568	i = bc_h(0)%i(m)
	569	j = bc_h(0)%j(m)
	570	k = bc_h(0)%k(m)
[2232]	571	tend(k-1,j,i) = tend(k,j,i)
[1]	572	ENDDO
[2232]	573	!
	574	!-- Downward facing
	575	!$OMP PARALLEL DO PRIVATE( i, j, k )
	576	DO m = 1, bc_h(1)%ns
[2298]	577	i = bc_h(1)%i(m)
	578	j = bc_h(1)%j(m)
	579	k = bc_h(1)%k(m)
[2232]	580	tend(k+1,j,i) = tend(k,j,i)
	581	ENDDO
[1]	582
	583	ELSE
	584	!
[2232]	585	!-- Dirichlet. Using surface data type, first for non-natural
	586	!-- surfaces, then for natural and urban surfaces
	587	!-- Upward facing
	588	!$OMP PARALLEL DO PRIVATE( i, j, k )
	589	DO m = 1, bc_h(0)%ns
[2298]	590	i = bc_h(0)%i(m)
	591	j = bc_h(0)%j(m)
	592	k = bc_h(0)%k(m)
[2232]	593	tend(k-1,j,i) = 0.0_wp
[1]	594	ENDDO
[2232]	595	!
	596	!-- Downward facing
	597	!$OMP PARALLEL DO PRIVATE( i, j, k )
	598	DO m = 1, bc_h(1)%ns
[2298]	599	i = bc_h(1)%i(m)
	600	j = bc_h(1)%j(m)
	601	k = bc_h(1)%k(m)
[2232]	602	tend(k+1,j,i) = 0.0_wp
	603	ENDDO
[1]	604
	605	ENDIF
	606
	607	!
	608	!-- Top boundary
	609	IF ( ibc_p_t == 1 ) THEN
	610	!
	611	!-- Neumann
[2232]	612	!$OMP PARALLEL DO PRIVATE (i,j,k)
[667]	613	DO i = nxlg, nxrg
	614	DO j = nysg, nyng
[1]	615	tend(nzt+1,j,i) = tend(nzt,j,i)
	616	ENDDO
	617	ENDDO
	618
	619	ELSE
	620	!
	621	!-- Dirichlet
[2232]	622	!$OMP PARALLEL DO PRIVATE (i,j,k)
[667]	623	DO i = nxlg, nxrg
	624	DO j = nysg, nyng
[1342]	625	tend(nzt+1,j,i) = 0.0_wp
[1]	626	ENDDO
	627	ENDDO
	628
	629	ENDIF
	630
	631	!
	632	!-- Exchange boundaries for p
[667]	633	CALL exchange_horiz( tend, nbgp )
[1]	634
	635	ELSEIF ( psolver == 'sor' ) THEN
	636
	637	!
	638	!-- Solve Poisson equation for perturbation pressure using SOR-Red/Black
	639	!-- scheme
[707]	640	CALL sor( d, ddzu_pres, ddzw, p_loc )
	641	tend = p_loc
[1]	642
[1575]	643	ELSEIF ( psolver(1:9) == 'multigrid' ) THEN
[1]	644
	645	!
	646	!-- Solve Poisson equation for perturbation pressure using Multigrid scheme,
[2298]	647	!-- array tend is used to store the residuals.
[778]	648
	649	!-- If the number of grid points of the gathered grid, which is collected
	650	!-- on PE0, is larger than the number of grid points of an PE, than array
	651	!-- tend will be enlarged.
	652	IF ( gathered_size > subdomain_size ) THEN
	653	DEALLOCATE( tend )
	654	ALLOCATE( tend(nzb:nzt_mg(mg_switch_to_pe0_level)+1,nys_mg( &
	655	mg_switch_to_pe0_level)-1:nyn_mg(mg_switch_to_pe0_level)+1,&
	656	nxl_mg(mg_switch_to_pe0_level)-1:nxr_mg( &
	657	mg_switch_to_pe0_level)+1) )
	658	ENDIF
	659
[1575]	660	IF ( psolver == 'multigrid' ) THEN
	661	CALL poismg( tend )
	662	ELSE
[1931]	663	CALL poismg_noopt( tend )
[1575]	664	ENDIF
[707]	665
[778]	666	IF ( gathered_size > subdomain_size ) THEN
	667	DEALLOCATE( tend )
	668	ALLOCATE( tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	669	ENDIF
	670
[1]	671	!
	672	!-- Restore perturbation pressure on tend because this array is used
	673	!-- further below to correct the velocity fields
[707]	674	DO i = nxl-1, nxr+1
	675	DO j = nys-1, nyn+1
	676	DO k = nzb, nzt+1
	677	tend(k,j,i) = p_loc(k,j,i)
	678	ENDDO
	679	ENDDO
	680	ENDDO
[667]	681
[1]	682	ENDIF
	683
	684	!
[707]	685	!-- Store perturbation pressure on array p, used for pressure data output.
	686	!-- Ghost layers are added in the output routines (except sor-method: see below)
[1918]	687	IF ( intermediate_timestep_count <= 1 ) THEN
[707]	688	!$OMP PARALLEL PRIVATE (i,j,k)
	689	!$OMP DO
	690	DO i = nxl-1, nxr+1
	691	DO j = nys-1, nyn+1
	692	DO k = nzb, nzt+1
	693	p(k,j,i) = tend(k,j,i) * &
[1929]	694	weight_substep_l
[673]	695	ENDDO
[1]	696	ENDDO
[707]	697	ENDDO
	698	!$OMP END PARALLEL
	699
[1762]	700	ELSEIF ( intermediate_timestep_count > 1 ) THEN
[707]	701	!$OMP PARALLEL PRIVATE (i,j,k)
	702	!$OMP DO
	703	DO i = nxl-1, nxr+1
	704	DO j = nys-1, nyn+1
	705	DO k = nzb, nzt+1
	706	p(k,j,i) = p(k,j,i) + tend(k,j,i) * &
[1929]	707	weight_substep_l
[673]	708	ENDDO
	709	ENDDO
[707]	710	ENDDO
	711	!$OMP END PARALLEL
	712
	713	ENDIF
[673]	714
[707]	715	!
	716	!-- SOR-method needs ghost layers for the next timestep
	717	IF ( psolver == 'sor' ) CALL exchange_horiz( p, nbgp )
[682]	718
[1]	719	!
	720	!-- Correction of the provisional velocities with the current perturbation
	721	!-- pressure just computed
[709]	722	IF ( conserve_volume_flow .AND. ( bc_lr_cyc .OR. bc_ns_cyc ) ) THEN
[1342]	723	volume_flow_l(1) = 0.0_wp
	724	volume_flow_l(2) = 0.0_wp
[1]	725	ENDIF
[3347]	726	!
	727	!-- Add pressure gradients to the velocity components. Note, the loops are
	728	!-- running over the entire model domain, even though, in case of non-cyclic
	729	!-- boundaries u- and v-component are not prognostic at i=0 and j=0,
	730	!-- respectiveley. However, in case of Dirichlet boundary conditions for the
	731	!-- velocities, zero-gradient conditions for the pressure are set, so that
	732	!-- no modification is imposed at the boundaries.
[1]	733	!$OMP PARALLEL PRIVATE (i,j,k)
	734	!$OMP DO
[673]	735	DO i = nxl, nxr
[1]	736	DO j = nys, nyn
[2118]	737
[2232]	738	DO k = nzb+1, nzt
	739	w(k,j,i) = w(k,j,i) - dt_3d * &
	740	( tend(k+1,j,i) - tend(k,j,i) ) * ddzu(k+1) &
	741	* weight_pres_l &
	742	* MERGE( 1.0_wp, 0.0_wp, &
	743	BTEST( wall_flags_0(k,j,i), 3 ) &
	744	)
[1]	745	ENDDO
[2118]	746
[2232]	747	DO k = nzb+1, nzt
	748	u(k,j,i) = u(k,j,i) - dt_3d * &
	749	( tend(k,j,i) - tend(k,j,i-1) ) * ddx &
	750	* weight_pres_l &
	751	* MERGE( 1.0_wp, 0.0_wp, &
	752	BTEST( wall_flags_0(k,j,i), 1 ) &
	753	)
[1]	754	ENDDO
[2118]	755
[2232]	756	DO k = nzb+1, nzt
	757	v(k,j,i) = v(k,j,i) - dt_3d * &
	758	( tend(k,j,i) - tend(k,j-1,i) ) * ddy &
	759	* weight_pres_l &
	760	* MERGE( 1.0_wp, 0.0_wp, &
	761	BTEST( wall_flags_0(k,j,i), 2 ) &
	762	)
[673]	763	ENDDO
[1]	764
	765	ENDDO
	766	ENDDO
	767	!$OMP END PARALLEL
[1113]	768
	769	!
[3012]	770	!-- The vertical velocity is not set to zero at nzt + 1 for nested domains
	771	!-- Instead it is set to the values of nzt (see routine vnest_boundary_conds
	772	!-- or pmci_interp_tril_t) BEFORE calling the pressure solver. To avoid jumps
	773	!-- while plotting profiles w at the top has to be set to the values in the
	774	!-- height nzt after above modifications. Hint: w level nzt+1 does not impact
	775	!-- results.
[3182]	776	IF ( child_domain .OR. coupling_mode == 'vnested_fine' ) THEN
[3012]	777	w(nzt+1,:,:) = w(nzt,:,:)
	778	ENDIF
	779
	780	!
[1113]	781	!-- Sum up the volume flow through the right and north boundary
[2232]	782	IF ( conserve_volume_flow .AND. bc_lr_cyc .AND. bc_ns_cyc .AND. &
[1113]	783	nxr == nx ) THEN
	784
	785	!$OMP PARALLEL PRIVATE (j,k)
	786	!$OMP DO
	787	DO j = nys, nyn
	788	!$OMP CRITICAL
[2232]	789	DO k = nzb+1, nzt
	790	volume_flow_l(1) = volume_flow_l(1) + u(k,j,nxr) * dzw(k) &
	791	* MERGE( 1.0_wp, 0.0_wp, &
	792	BTEST( wall_flags_0(k,j,nxr), 1 )&
	793	)
[1113]	794	ENDDO
	795	!$OMP END CRITICAL
	796	ENDDO
	797	!$OMP END PARALLEL
	798
	799	ENDIF
	800
[2232]	801	IF ( conserve_volume_flow .AND. bc_ns_cyc .AND. bc_lr_cyc .AND. &
[1113]	802	nyn == ny ) THEN
	803
	804	!$OMP PARALLEL PRIVATE (i,k)
	805	!$OMP DO
	806	DO i = nxl, nxr
	807	!$OMP CRITICAL
[2232]	808	DO k = nzb+1, nzt
	809	volume_flow_l(2) = volume_flow_l(2) + v(k,nyn,i) * dzw(k) &
	810	* MERGE( 1.0_wp, 0.0_wp, &
	811	BTEST( wall_flags_0(k,nyn,i), 2 )&
	812	)
[1113]	813	ENDDO
	814	!$OMP END CRITICAL
	815	ENDDO
	816	!$OMP END PARALLEL
	817
	818	ENDIF
[673]	819
[1]	820	!
	821	!-- Conserve the volume flow
[707]	822	IF ( conserve_volume_flow .AND. ( bc_lr_cyc .AND. bc_ns_cyc ) ) THEN
[1]	823
	824	#if defined( __parallel )
[622]	825	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[1]	826	CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 2, MPI_REAL, &
	827	MPI_SUM, comm2d, ierr )
	828	#else
	829	volume_flow = volume_flow_l
	830	#endif
	831
[1799]	832	volume_flow_offset(1:2) = ( volume_flow_initial(1:2) - volume_flow(1:2) ) / &
	833	volume_flow_area(1:2)
[1]	834
	835	!$OMP PARALLEL PRIVATE (i,j,k)
	836	!$OMP DO
	837	DO i = nxl, nxr
	838	DO j = nys, nyn
[2232]	839	DO k = nzb+1, nzt
	840	u(k,j,i) = u(k,j,i) + volume_flow_offset(1) &
	841	* MERGE( 1.0_wp, 0.0_wp, &
	842	BTEST( wall_flags_0(k,j,i), 1 ) &
	843	)
[719]	844	ENDDO
[2232]	845	DO k = nzb+1, nzt
	846	v(k,j,i) = v(k,j,i) + volume_flow_offset(2) &
	847	* MERGE( 1.0_wp, 0.0_wp, &
	848	BTEST( wall_flags_0(k,j,i), 2 ) &
	849	)
[667]	850	ENDDO
[1]	851	ENDDO
	852	ENDDO
[667]	853
[1]	854	!$OMP END PARALLEL
	855
	856	ENDIF
	857
	858	!
	859	!-- Exchange of boundaries for the velocities
[667]	860	CALL exchange_horiz( u, nbgp )
	861	CALL exchange_horiz( v, nbgp )
	862	CALL exchange_horiz( w, nbgp )
[1]	863
	864	!
	865	!-- Compute the divergence of the corrected velocity field,
[1908]	866	!-- A possible PE-sum is computed in flow_statistics. Carry out computation
	867	!-- only at last Runge-Kutta step.
[1918]	868	IF ( intermediate_timestep_count == intermediate_timestep_count_max .OR. &
	869	intermediate_timestep_count == 0 ) THEN
[1908]	870	CALL cpu_log( log_point_s(1), 'divergence', 'start' )
	871	sums_divnew_l = 0.0_wp
[1]	872
	873	!
[1908]	874	!-- d must be reset to zero because it can contain nonzero values below the
	875	!-- topography
	876	IF ( topography /= 'flat' ) d = 0.0_wp
[1]	877
[1908]	878	localsum = 0.0_wp
	879	threadsum = 0.0_wp
[1]	880
[1908]	881	!$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
[2073]	882	#if defined( __ibm )
[1908]	883	!$OMP DO SCHEDULE( STATIC )
	884	DO i = nxl, nxr
	885	DO j = nys, nyn
[2232]	886	DO k = nzb+1, nzt
	887	d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx + &
	888	( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy + &
	889	( w(k,j,i) * rho_air_zw(k) - &
	890	w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k) &
	891	) * MERGE( 1.0_wp, 0.0_wp, &
	892	BTEST( wall_flags_0(k,j,i), 0 ) &
	893	)
[1908]	894	ENDDO
	895	DO k = nzb+1, nzt
	896	threadsum = threadsum + ABS( d(k,j,i) )
	897	ENDDO
[1]	898	ENDDO
	899	ENDDO
	900	#else
[2073]	901	!$OMP DO SCHEDULE( STATIC )
[1908]	902	DO i = nxl, nxr
	903	DO j = nys, nyn
[2232]	904	DO k = nzb+1, nzt
[2037]	905	d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx + &
	906	( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy + &
	907	( w(k,j,i) * rho_air_zw(k) - &
	908	w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k) &
[2232]	909	) * MERGE( 1.0_wp, 0.0_wp, &
	910	BTEST( wall_flags_0(k,j,i), 0 ) &
	911	)
[1908]	912	ENDDO
[1113]	913	ENDDO
	914	ENDDO
	915	!
[1908]	916	!-- Compute possible PE-sum of divergences for flow_statistics
[2073]	917	!$OMP DO SCHEDULE( STATIC )
[1908]	918	DO i = nxl, nxr
	919	DO j = nys, nyn
	920	DO k = nzb+1, nzt
	921	threadsum = threadsum + ABS( d(k,j,i) )
	922	ENDDO
[1]	923	ENDDO
	924	ENDDO
	925	#endif
[667]	926
[1908]	927	localsum = localsum + threadsum
	928	!$OMP END PARALLEL
[1]	929
	930	!
[1908]	931	!-- For completeness, set the divergence sum of all statistic regions to those
	932	!-- of the total domain
	933	sums_divnew_l(0:statistic_regions) = localsum
[1]	934
[1908]	935	CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
[1]	936
[1908]	937	ENDIF
	938
[1]	939	CALL cpu_log( log_point(8), 'pres', 'stop' )
	940
	941
	942	END SUBROUTINE pres

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |