Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |