Home

Context Navigation

source: palm/trunk/SOURCE/boundary_conds.f90 @ 4114

Last change on this file since 4114 was 4102, checked in by suehring, 5 years ago
Bugfix, set Neumann boundary conditions for the subgrid TKE at vertical walls instead of implicit Dirichlet conditions that always act as a sink term for the subgrid TKE. Therefore, add new data structure for vertical surfaces and revise the setting of the boundary grid point index space. Moreover, slightly revise setting of boundary conditions at upward- and downward facing surfaces. Finally, setting of boundary conditions for subgrid TKE and dissipation (in RANS mode) is now modularized. Update test case results.
Property svn:keywords set to `Id`
File size: 40.6 KB

Rev	Line
[1682]	1	!> @file boundary_conds.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	!
[3655]	17	! Copyright 1997-2019 Leibniz Universitaet Hannover
[2000]	18	!------------------------------------------------------------------------------!
[1036]	19	!
[484]	20	! Current revisions:
[1]	21	! -----------------
[1933]	22	!
[3589]	23	!
[1321]	24	! Former revisions:
	25	! -----------------
	26	! $Id: boundary_conds.f90 4102 2019-07-17 16:00:03Z schwenkel $
[4102]	27	! - Revise setting for boundary conditions at horizontal walls, use the offset
	28	! index that belongs to the data structure instead of pre-calculating
	29	! the offset index for each facing.
	30	! - Set boundary conditions for bulk-cloud quantities also at downward facing
	31	! surfaces
	32	!
	33	! 4087 2019-07-11 11:35:23Z gronemeier
[4087]	34	! Add comment
	35	!
	36	! 4086 2019-07-11 05:55:44Z gronemeier
[4086]	37	! Bugfix: use constant-flux layer condition for e in all rans modes
	38	!
	39	! 3879 2019-04-08 20:25:23Z knoop
[3717]	40	! Bugfix, do not set boundary conditions for potential temperature in neutral
	41	! runs.
	42	!
	43	! 3655 2019-01-07 16:51:22Z knoop
[3634]	44	! OpenACC port for SPEC
	45	!
	46	! 3589 2018-11-30 15:09:51Z suehring
[3589]	47	! Move the control parameter "salsa" from salsa_mod to control_parameters
	48	! (M. Kurppa)
	49	!
	50	! 3582 2018-11-29 19:16:36Z suehring
[3562]	51	! variables documented
	52	!
	53	! 3467 2018-10-30 19:05:21Z suehring
[3467]	54	! Implementation of a new aerosol module salsa.
	55	!
	56	! 3347 2018-10-15 14:21:08Z suehring
[3347]	57	! Bugfix in setting boundary conditions in offline nesting
	58	!
	59	! 3341 2018-10-15 10:31:27Z suehring
[3294]	60	! changes concerning modularization of ocean option
	61	!
	62	! 3274 2018-09-24 15:42:55Z knoop
[3274]	63	! Modularization of all bulk cloud physics code components
	64	!
	65	! 3241 2018-09-12 15:02:00Z raasch
[3241]	66	! unused variables removed
	67	!
	68	! 3183 2018-07-27 14:25:55Z suehring
[3183]	69	! Rename some variables concerning LES-LES as well as offline nesting
	70	!
	71	! 3182 2018-07-27 13:36:03Z suehring
[3129]	72	! - Use wall function for e_p and diss_p in case of rans_tke_e
	73	! - move limitation of diss_p from tcm_prognostic
	74	!
	75	! 2938 2018-03-27 15:52:42Z suehring
[2938]	76	! Set boundary condition for TKE and TKE dissipation rate in case of nesting
	77	! and if parent model operates in RANS mode but child model in LES mode.
	78	! mode
	79	!
	80	! 2793 2018-02-07 10:54:33Z suehring
[2766]	81	! Removed preprocessor directive __chem
	82	!
	83	! 2718 2018-01-02 08:49:38Z maronga
[2716]	84	! Corrected "Former revisions" section
	85	!
	86	! 2696 2017-12-14 17:12:51Z kanani
	87	! Change in file header (GPL part)
[2696]	88	! Adjust boundary conditions for e and diss in case of TKE-e closure (TG)
	89	! Implementation of chemistry module (FK)
	90	!
	91	! 2569 2017-10-20 11:54:42Z kanani
[2569]	92	! Removed redundant code for ibc_s_b=1 and ibc_q_b=1
	93	!
	94	! 2365 2017-08-21 14:59:59Z kanani
[2365]	95	! Vertical grid nesting implemented: exclude setting vertical velocity to zero
	96	! on fine grid (SadiqHuq)
	97	!
	98	! 2320 2017-07-21 12:47:43Z suehring
[2320]	99	! Remove unused control parameter large_scale_forcing from only-list
	100	!
	101	! 2292 2017-06-20 09:51:42Z schwenkel
[2292]	102	! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
	103	! includes two more prognostic equations for cloud drop concentration (nc)
	104	! and cloud water content (qc).
	105	!
	106	! 2233 2017-05-30 18:08:54Z suehring
[1321]	107	!
[2233]	108	! 2232 2017-05-30 17:47:52Z suehring
	109	! Set boundary conditions on topography top using flag method.
	110	!
[2119]	111	! 2118 2017-01-17 16:38:49Z raasch
	112	! OpenACC directives removed
	113	!
[2001]	114	! 2000 2016-08-20 18:09:15Z knoop
	115	! Forced header and separation lines into 80 columns
	116	!
[1993]	117	! 1992 2016-08-12 15:14:59Z suehring
	118	! Adjustments for top boundary condition for passive scalar
	119	!
[1961]	120	! 1960 2016-07-12 16:34:24Z suehring
	121	! Treat humidity and passive scalar separately
	122	!
[1933]	123	! 1823 2016-04-07 08:57:52Z hoffmann
	124	! Initial version of purely vertical nesting introduced.
	125	!
[1823]	126	! 1822 2016-04-07 07:49:42Z hoffmann
	127	! icloud_scheme removed. microphyisics_seifert added.
	128	!
[1765]	129	! 1764 2016-02-28 12:45:19Z raasch
	130	! index bug for u_p at left outflow removed
	131	!
[1763]	132	! 1762 2016-02-25 12:31:13Z hellstea
	133	! Introduction of nested domain feature
	134	!
[1744]	135	! 1742 2016-01-13 09:50:06Z raasch
	136	! bugfix for outflow Neumann boundary conditions at bottom and top
	137	!
[1718]	138	! 1717 2015-11-11 15:09:47Z raasch
	139	! Bugfix: index error in outflow conditions for left boundary
	140	!
[1683]	141	! 1682 2015-10-07 23:56:08Z knoop
	142	! Code annotations made doxygen readable
	143	!
[1717]	144	! 1410 2014-05-23 12:16:18Z suehring
[1463]	145	! Bugfix: set dirichlet boundary condition for passive_scalar at model domain
	146	! top
	147	!
[1410]	148	! 1399 2014-05-07 11:16:25Z heinze
	149	! Bugfix: set inflow boundary conditions also if no humidity or passive_scalar
	150	! is used.
	151	!
[1399]	152	! 1398 2014-05-07 11:15:00Z heinze
	153	! Dirichlet-condition at the top for u and v changed to u_init and v_init also
	154	! for large_scale_forcing
	155	!
[1381]	156	! 1380 2014-04-28 12:40:45Z heinze
	157	! Adjust Dirichlet-condition at the top for pt in case of nudging
	158	!
[1362]	159	! 1361 2014-04-16 15:17:48Z hoffmann
	160	! Bottom and top boundary conditions of rain water content (qr) and
	161	! rain drop concentration (nr) changed to Dirichlet
	162	!
[1354]	163	! 1353 2014-04-08 15:21:23Z heinze
	164	! REAL constants provided with KIND-attribute
	165	!
[1321]	166	! 1320 2014-03-20 08:40:49Z raasch
[1320]	167	! ONLY-attribute added to USE-statements,
	168	! kind-parameters added to all INTEGER and REAL declaration statements,
	169	! kinds are defined in new module kinds,
	170	! revision history before 2012 removed,
	171	! comment fields (!:) to be used for variable explanations added to
	172	! all variable declaration statements
[1160]	173	!
[1258]	174	! 1257 2013-11-08 15:18:40Z raasch
	175	! loop independent clauses added
	176	!
[1242]	177	! 1241 2013-10-30 11:36:58Z heinze
	178	! Adjust ug and vg at each timestep in case of large_scale_forcing
	179	!
[1160]	180	! 1159 2013-05-21 11:58:22Z fricke
[1159]	181	! Bugfix: Neumann boundary conditions for the velocity components at the
	182	! outflow are in fact radiation boundary conditions using the maximum phase
	183	! velocity that ensures numerical stability (CFL-condition).
	184	! Hence, logical operator use_cmax is now used instead of bc_lr_dirneu/_neudir.
	185	! Bugfix: In case of use_cmax at the outflow, u, v, w are replaced by
	186	! u_p, v_p, w_p
[1116]	187	!
	188	! 1115 2013-03-26 18:16:16Z hoffmann
	189	! boundary conditions of two-moment cloud scheme are restricted to Neumann-
	190	! boundary-conditions
	191	!
[1114]	192	! 1113 2013-03-10 02:48:14Z raasch
	193	! GPU-porting
	194	! dummy argument "range" removed
	195	! Bugfix: wrong index in loops of radiation boundary condition
[1113]	196	!
[1054]	197	! 1053 2012-11-13 17:11:03Z hoffmann
	198	! boundary conditions for the two new prognostic equations (nr, qr) of the
	199	! two-moment cloud scheme
	200	!
[1037]	201	! 1036 2012-10-22 13:43:42Z raasch
	202	! code put under GPL (PALM 3.9)
	203	!
[997]	204	! 996 2012-09-07 10:41:47Z raasch
	205	! little reformatting
	206	!
[979]	207	! 978 2012-08-09 08:28:32Z fricke
	208	! Neumann boudnary conditions are added at the inflow boundary for the SGS-TKE.
	209	! Outflow boundary conditions for the velocity components can be set to Neumann
	210	! conditions or to radiation conditions with a horizontal averaged phase
	211	! velocity.
	212	!
[876]	213	! 875 2012-04-02 15:35:15Z gryschka
	214	! Bugfix in case of dirichlet inflow bc at the right or north boundary
	215	!
[1]	216	! Revision 1.1 1997/09/12 06:21:34 raasch
	217	! Initial revision
	218	!
	219	!
	220	! Description:
	221	! ------------
[1682]	222	!> Boundary conditions for the prognostic quantities.
	223	!> One additional bottom boundary condition is applied for the TKE (=(u)*2)
	224	!> in prandtl_fluxes. The cyclic lateral boundary conditions are implicitly
	225	!> handled in routine exchange_horiz. Pressure boundary conditions are
	226	!> explicitly set in routines pres, poisfft, poismg and sor.
[1]	227	!------------------------------------------------------------------------------!
[1682]	228	SUBROUTINE boundary_conds
	229
[1]	230
[1320]	231	USE arrays_3d, &
	232	ONLY: c_u, c_u_m, c_u_m_l, c_v, c_v_m, c_v_m_l, c_w, c_w_m, c_w_m_l, &
[4102]	233	dzu, nc_p, nr_p, pt, pt_init, pt_p, q, &
[3241]	234	q_p, qc_p, qr_p, s, s_p, sa, sa_p, u, u_init, u_m_l, u_m_n, &
	235	u_m_r, u_m_s, u_p, v, v_init, v_m_l, v_m_n, v_m_r, v_m_s, v_p, &
	236	w, w_p, w_m_l, w_m_n, w_m_r, w_m_s
[2696]	237
[3294]	238	USE bulk_cloud_model_mod, &
	239	ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert
	240
[2696]	241	USE chemistry_model_mod, &
[3879]	242	ONLY: chem_boundary_conds
	243
[1320]	244	USE control_parameters, &
[4102]	245	ONLY: air_chemistry, bc_dirichlet_l, &
[3182]	246	bc_dirichlet_s, bc_radiation_l, bc_radiation_n, bc_radiation_r, &
	247	bc_radiation_s, bc_pt_t_val, bc_q_t_val, bc_s_t_val, &
[4102]	248	child_domain, coupling_mode, dt_3d, &
[3582]	249	humidity, ibc_pt_b, ibc_pt_t, ibc_q_b, ibc_q_t, ibc_s_b, &
	250	ibc_s_t, ibc_uv_b, ibc_uv_t, intermediate_timestep_count, &
[3717]	251	nesting_offline, neutral, nudging, ocean_mode, passive_scalar, &
[4102]	252	tsc, salsa, use_cmax
[1320]	253
	254	USE grid_variables, &
	255	ONLY: ddx, ddy, dx, dy
	256
	257	USE indices, &
[3294]	258	ONLY: nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb, nzt
[1320]	259
	260	USE kinds
	261
[3294]	262	USE ocean_mod, &
	263	ONLY: ibc_sa_t
[3274]	264
[1]	265	USE pegrid
	266
[1933]	267	USE pmc_interface, &
[4102]	268	ONLY : nesting_mode
[3467]	269
	270	USE salsa_mod, &
[3582]	271	ONLY: salsa_boundary_conds
[1320]	272
[2232]	273	USE surface_mod, &
[4102]	274	ONLY : bc_h
[1933]	275
[3129]	276	USE turbulence_closure_mod, &
[4102]	277	ONLY: tcm_boundary_conds
[3129]	278
[1]	279	IMPLICIT NONE
	280
[2232]	281	INTEGER(iwp) :: i !< grid index x direction
	282	INTEGER(iwp) :: j !< grid index y direction
	283	INTEGER(iwp) :: k !< grid index z direction
	284	INTEGER(iwp) :: l !< running index boundary type, for up- and downward-facing walls
	285	INTEGER(iwp) :: m !< running index surface elements
[1]	286
[3562]	287	REAL(wp) :: c_max !< maximum phase velocity allowed by CFL criterion, used for outflow boundary condition
	288	REAL(wp) :: denom !< horizontal gradient of velocity component normal to the outflow boundary
[1]	289
[73]	290
[1]	291	!
[1113]	292	!-- Bottom boundary
	293	IF ( ibc_uv_b == 1 ) THEN
	294	u_p(nzb,:,:) = u_p(nzb+1,:,:)
	295	v_p(nzb,:,:) = v_p(nzb+1,:,:)
	296	ENDIF
[2232]	297	!
	298	!-- Set zero vertical velocity at topography top (l=0), or bottom (l=1) in case
	299	!-- of downward-facing surfaces.
	300	DO l = 0, 1
	301	!$OMP PARALLEL DO PRIVATE( i, j, k )
[3634]	302	!$ACC PARALLEL LOOP PRIVATE(i, j, k) &
	303	!$ACC PRESENT(bc_h, w_p)
[2232]	304	DO m = 1, bc_h(l)%ns
	305	i = bc_h(l)%i(m)
	306	j = bc_h(l)%j(m)
	307	k = bc_h(l)%k(m)
[4102]	308	w_p(k+bc_h(l)%koff,j,i) = 0.0_wp
[1113]	309	ENDDO
	310	ENDDO
	311
	312	!
[1762]	313	!-- Top boundary. A nested domain ( ibc_uv_t = 3 ) does not require settings.
[1113]	314	IF ( ibc_uv_t == 0 ) THEN
[3634]	315	!$ACC KERNELS PRESENT(u_p, v_p, u_init, v_init)
[1113]	316	u_p(nzt+1,:,:) = u_init(nzt+1)
	317	v_p(nzt+1,:,:) = v_init(nzt+1)
[3634]	318	!$ACC END KERNELS
[1762]	319	ELSEIF ( ibc_uv_t == 1 ) THEN
[1113]	320	u_p(nzt+1,:,:) = u_p(nzt,:,:)
	321	v_p(nzt+1,:,:) = v_p(nzt,:,:)
	322	ENDIF
	323
[2365]	324	!
	325	!-- Vertical nesting: Vertical velocity not zero at the top of the fine grid
[3347]	326	IF ( .NOT. child_domain .AND. .NOT. nesting_offline .AND. &
[2365]	327	TRIM(coupling_mode) /= 'vnested_fine' ) THEN
[3634]	328	!$ACC KERNELS PRESENT(w_p)
[2365]	329	w_p(nzt:nzt+1,:,:) = 0.0_wp !< nzt is not a prognostic level (but cf. pres)
[3634]	330	!$ACC END KERNELS
[1762]	331	ENDIF
	332
[1113]	333	!
[2232]	334	!-- Temperature at bottom and top boundary.
[1113]	335	!-- In case of coupled runs (ibc_pt_b = 2) the temperature is given by
	336	!-- the sea surface temperature of the coupled ocean model.
[2232]	337	!-- Dirichlet
[3717]	338	IF ( .NOT. neutral ) THEN
	339	IF ( ibc_pt_b == 0 ) THEN
	340	DO l = 0, 1
	341	!$OMP PARALLEL DO PRIVATE( i, j, k )
	342	DO m = 1, bc_h(l)%ns
	343	i = bc_h(l)%i(m)
	344	j = bc_h(l)%j(m)
	345	k = bc_h(l)%k(m)
[4102]	346	pt_p(k+bc_h(l)%koff,j,i) = pt(k+bc_h(l)%koff,j,i)
[3717]	347	ENDDO
[1]	348	ENDDO
[3717]	349	!
	350	!-- Neumann, zero-gradient
	351	ELSEIF ( ibc_pt_b == 1 ) THEN
	352	DO l = 0, 1
	353	!$OMP PARALLEL DO PRIVATE( i, j, k )
	354	!$ACC PARALLEL LOOP PRIVATE(i, j, k) &
	355	!$ACC PRESENT(bc_h, pt_p)
	356	DO m = 1, bc_h(l)%ns
	357	i = bc_h(l)%i(m)
	358	j = bc_h(l)%j(m)
	359	k = bc_h(l)%k(m)
[4102]	360	pt_p(k+bc_h(l)%koff,j,i) = pt_p(k,j,i)
[3717]	361	ENDDO
[1113]	362	ENDDO
[3717]	363	ENDIF
	364
	365	!
	366	!-- Temperature at top boundary
	367	IF ( ibc_pt_t == 0 ) THEN
	368	pt_p(nzt+1,:,:) = pt(nzt+1,:,:)
	369	!
	370	!-- In case of nudging adjust top boundary to pt which is
	371	!-- read in from NUDGING-DATA
	372	IF ( nudging ) THEN
	373	pt_p(nzt+1,:,:) = pt_init(nzt+1)
	374	ENDIF
	375	ELSEIF ( ibc_pt_t == 1 ) THEN
	376	pt_p(nzt+1,:,:) = pt_p(nzt,:,:)
	377	ELSEIF ( ibc_pt_t == 2 ) THEN
	378	!$ACC KERNELS PRESENT(pt_p, dzu)
	379	pt_p(nzt+1,:,:) = pt_p(nzt,:,:) + bc_pt_t_val * dzu(nzt+1)
	380	!$ACC END KERNELS
	381	ENDIF
[1113]	382	ENDIF
[1]	383	!
[1113]	384	!-- Boundary conditions for salinity
[3294]	385	IF ( ocean_mode ) THEN
[1113]	386	!
	387	!-- Bottom boundary: Neumann condition because salinity flux is always
[2232]	388	!-- given.
	389	DO l = 0, 1
	390	!$OMP PARALLEL DO PRIVATE( i, j, k )
	391	DO m = 1, bc_h(l)%ns
	392	i = bc_h(l)%i(m)
	393	j = bc_h(l)%j(m)
	394	k = bc_h(l)%k(m)
[4102]	395	sa_p(k+bc_h(l)%koff,j,i) = sa_p(k,j,i)
[1]	396	ENDDO
[1113]	397	ENDDO
[1]	398	!
[1113]	399	!-- Top boundary: Dirichlet or Neumann
	400	IF ( ibc_sa_t == 0 ) THEN
	401	sa_p(nzt+1,:,:) = sa(nzt+1,:,:)
	402	ELSEIF ( ibc_sa_t == 1 ) THEN
	403	sa_p(nzt+1,:,:) = sa_p(nzt,:,:)
[1]	404	ENDIF
	405
[1113]	406	ENDIF
	407
[1]	408	!
[1960]	409	!-- Boundary conditions for total water content,
[1113]	410	!-- bottom and top boundary (see also temperature)
[1960]	411	IF ( humidity ) THEN
[1113]	412	!
	413	!-- Surface conditions for constant_humidity_flux
[2232]	414	!-- Run loop over all non-natural and natural walls. Note, in wall-datatype
	415	!-- the k coordinate belongs to the atmospheric grid point, therefore, set
	416	!-- q_p at k-1
[1113]	417	IF ( ibc_q_b == 0 ) THEN
[2232]	418
	419	DO l = 0, 1
	420	!$OMP PARALLEL DO PRIVATE( i, j, k )
	421	DO m = 1, bc_h(l)%ns
	422	i = bc_h(l)%i(m)
	423	j = bc_h(l)%j(m)
	424	k = bc_h(l)%k(m)
[4102]	425	q_p(k+bc_h(l)%koff,j,i) = q(k+bc_h(l)%koff,j,i)
[1]	426	ENDDO
	427	ENDDO
[2232]	428
[1113]	429	ELSE
[2232]	430
	431	DO l = 0, 1
	432	!$OMP PARALLEL DO PRIVATE( i, j, k )
	433	DO m = 1, bc_h(l)%ns
	434	i = bc_h(l)%i(m)
	435	j = bc_h(l)%j(m)
	436	k = bc_h(l)%k(m)
[4102]	437	q_p(k+bc_h(l)%koff,j,i) = q_p(k,j,i)
[95]	438	ENDDO
	439	ENDDO
[1113]	440	ENDIF
[95]	441	!
[1113]	442	!-- Top boundary
[1462]	443	IF ( ibc_q_t == 0 ) THEN
	444	q_p(nzt+1,:,:) = q(nzt+1,:,:)
	445	ELSEIF ( ibc_q_t == 1 ) THEN
[1992]	446	q_p(nzt+1,:,:) = q_p(nzt,:,:) + bc_q_t_val * dzu(nzt+1)
[1462]	447	ENDIF
[95]	448
[3274]	449	IF ( bulk_cloud_model .AND. microphysics_morrison ) THEN
[2292]	450	!
	451	!-- Surface conditions cloud water (Dirichlet)
	452	!-- Run loop over all non-natural and natural walls. Note, in wall-datatype
	453	!-- the k coordinate belongs to the atmospheric grid point, therefore, set
[4102]	454	!-- qr_p and nr_p at upward (k-1) and downward-facing (k+1) walls
	455	DO l = 0, 1
[2292]	456	!$OMP PARALLEL DO PRIVATE( i, j, k )
[4102]	457	DO m = 1, bc_h(l)%ns
	458	i = bc_h(l)%i(m)
	459	j = bc_h(l)%j(m)
	460	k = bc_h(l)%k(m)
	461	qc_p(k+bc_h(l)%koff,j,i) = 0.0_wp
	462	nc_p(k+bc_h(l)%koff,j,i) = 0.0_wp
	463	ENDDO
[2292]	464	ENDDO
	465	!
	466	!-- Top boundary condition for cloud water (Dirichlet)
	467	qc_p(nzt+1,:,:) = 0.0_wp
	468	nc_p(nzt+1,:,:) = 0.0_wp
	469
	470	ENDIF
	471
[3274]	472	IF ( bulk_cloud_model .AND. microphysics_seifert ) THEN
[1113]	473	!
[1361]	474	!-- Surface conditions rain water (Dirichlet)
[2232]	475	!-- Run loop over all non-natural and natural walls. Note, in wall-datatype
	476	!-- the k coordinate belongs to the atmospheric grid point, therefore, set
[4102]	477	!-- qr_p and nr_p at upward (k-1) and downward-facing (k+1) walls
	478	DO l = 0, 1
[2232]	479	!$OMP PARALLEL DO PRIVATE( i, j, k )
[4102]	480	DO m = 1, bc_h(l)%ns
	481	i = bc_h(l)%i(m)
	482	j = bc_h(l)%j(m)
	483	k = bc_h(l)%k(m)
	484	qr_p(k+bc_h(l)%koff,j,i) = 0.0_wp
	485	nr_p(k+bc_h(l)%koff,j,i) = 0.0_wp
	486	ENDDO
[1115]	487	ENDDO
[1]	488	!
[1361]	489	!-- Top boundary condition for rain water (Dirichlet)
	490	qr_p(nzt+1,:,:) = 0.0_wp
	491	nr_p(nzt+1,:,:) = 0.0_wp
[1115]	492
[1]	493	ENDIF
[1409]	494	ENDIF
[1]	495	!
[1960]	496	!-- Boundary conditions for scalar,
	497	!-- bottom and top boundary (see also temperature)
	498	IF ( passive_scalar ) THEN
	499	!
	500	!-- Surface conditions for constant_humidity_flux
[2232]	501	!-- Run loop over all non-natural and natural walls. Note, in wall-datatype
	502	!-- the k coordinate belongs to the atmospheric grid point, therefore, set
	503	!-- s_p at k-1
[1960]	504	IF ( ibc_s_b == 0 ) THEN
[2232]	505
	506	DO l = 0, 1
	507	!$OMP PARALLEL DO PRIVATE( i, j, k )
	508	DO m = 1, bc_h(l)%ns
	509	i = bc_h(l)%i(m)
	510	j = bc_h(l)%j(m)
	511	k = bc_h(l)%k(m)
[4102]	512	s_p(k+bc_h(l)%koff,j,i) = s(k+bc_h(l)%koff,j,i)
[1960]	513	ENDDO
	514	ENDDO
[2232]	515
[1960]	516	ELSE
[2232]	517
	518	DO l = 0, 1
	519	!$OMP PARALLEL DO PRIVATE( i, j, k )
	520	DO m = 1, bc_h(l)%ns
	521	i = bc_h(l)%i(m)
	522	j = bc_h(l)%j(m)
	523	k = bc_h(l)%k(m)
[4102]	524	s_p(k+bc_h(l)%koff,j,i) = s_p(k,j,i)
[1960]	525	ENDDO
	526	ENDDO
	527	ENDIF
	528	!
[1992]	529	!-- Top boundary condition
	530	IF ( ibc_s_t == 0 ) THEN
[1960]	531	s_p(nzt+1,:,:) = s(nzt+1,:,:)
[1992]	532	ELSEIF ( ibc_s_t == 1 ) THEN
	533	s_p(nzt+1,:,:) = s_p(nzt,:,:)
	534	ELSEIF ( ibc_s_t == 2 ) THEN
	535	s_p(nzt+1,:,:) = s_p(nzt,:,:) + bc_s_t_val * dzu(nzt+1)
[1960]	536	ENDIF
	537
[4102]	538	ENDIF
[1960]	539	!
[4102]	540	!-- Set boundary conditions for subgrid TKE and dissipation (RANS mode)
	541	CALL tcm_boundary_conds
	542	!
[2696]	543	!-- Top/bottom boundary conditions for chemical species
	544	IF ( air_chemistry ) CALL chem_boundary_conds( 'set_bc_bottomtop' )
	545	!
[1762]	546	!-- In case of inflow or nest boundary at the south boundary the boundary for v
	547	!-- is at nys and in case of inflow or nest boundary at the left boundary the
	548	!-- boundary for u is at nxl. Since in prognostic_equations (cache optimized
	549	!-- version) these levels are handled as a prognostic level, boundary values
	550	!-- have to be restored here.
[3182]	551	IF ( bc_dirichlet_s ) THEN
[1409]	552	v_p(:,nys,:) = v_p(:,nys-1,:)
[3182]	553	ELSEIF ( bc_dirichlet_l ) THEN
[1409]	554	u_p(:,:,nxl) = u_p(:,:,nxl-1)
	555	ENDIF
[1]	556
	557	!
[1762]	558	!-- The same restoration for u at i=nxl and v at j=nys as above must be made
[1933]	559	!-- in case of nest boundaries. This must not be done in case of vertical nesting
[3182]	560	!-- mode as in that case the lateral boundaries are actually cyclic.
[4102]	561	!-- Lateral oundary conditions for TKE and dissipation are set
	562	!-- in tcm_boundary_conds.
[3182]	563	IF ( nesting_mode /= 'vertical' .OR. nesting_offline ) THEN
	564	IF ( bc_dirichlet_s ) THEN
[1933]	565	v_p(:,nys,:) = v_p(:,nys-1,:)
	566	ENDIF
[3182]	567	IF ( bc_dirichlet_l ) THEN
[1933]	568	u_p(:,:,nxl) = u_p(:,:,nxl-1)
	569	ENDIF
[1762]	570	ENDIF
	571
	572	!
[4102]	573	!-- Lateral boundary conditions for scalar quantities at the outflow.
	574	!-- Lateral oundary conditions for TKE and dissipation are set
	575	!-- in tcm_boundary_conds.
[3182]	576	IF ( bc_radiation_s ) THEN
[1409]	577	pt_p(:,nys-1,:) = pt_p(:,nys,:)
[1960]	578	IF ( humidity ) THEN
[1409]	579	q_p(:,nys-1,:) = q_p(:,nys,:)
[3274]	580	IF ( bulk_cloud_model .AND. microphysics_morrison ) THEN
[2292]	581	qc_p(:,nys-1,:) = qc_p(:,nys,:)
	582	nc_p(:,nys-1,:) = nc_p(:,nys,:)
	583	ENDIF
[3274]	584	IF ( bulk_cloud_model .AND. microphysics_seifert ) THEN
[1409]	585	qr_p(:,nys-1,:) = qr_p(:,nys,:)
	586	nr_p(:,nys-1,:) = nr_p(:,nys,:)
[1053]	587	ENDIF
[1409]	588	ENDIF
[1960]	589	IF ( passive_scalar ) s_p(:,nys-1,:) = s_p(:,nys,:)
[3182]	590	ELSEIF ( bc_radiation_n ) THEN
[1409]	591	pt_p(:,nyn+1,:) = pt_p(:,nyn,:)
[1960]	592	IF ( humidity ) THEN
[1409]	593	q_p(:,nyn+1,:) = q_p(:,nyn,:)
[3274]	594	IF ( bulk_cloud_model .AND. microphysics_morrison ) THEN
[2292]	595	qc_p(:,nyn+1,:) = qc_p(:,nyn,:)
	596	nc_p(:,nyn+1,:) = nc_p(:,nyn,:)
	597	ENDIF
[3274]	598	IF ( bulk_cloud_model .AND. microphysics_seifert ) THEN
[1409]	599	qr_p(:,nyn+1,:) = qr_p(:,nyn,:)
	600	nr_p(:,nyn+1,:) = nr_p(:,nyn,:)
[1053]	601	ENDIF
[1409]	602	ENDIF
[1960]	603	IF ( passive_scalar ) s_p(:,nyn+1,:) = s_p(:,nyn,:)
[3182]	604	ELSEIF ( bc_radiation_l ) THEN
[1409]	605	pt_p(:,:,nxl-1) = pt_p(:,:,nxl)
[1960]	606	IF ( humidity ) THEN
[1409]	607	q_p(:,:,nxl-1) = q_p(:,:,nxl)
[3274]	608	IF ( bulk_cloud_model .AND. microphysics_morrison ) THEN
[2292]	609	qc_p(:,:,nxl-1) = qc_p(:,:,nxl)
	610	nc_p(:,:,nxl-1) = nc_p(:,:,nxl)
	611	ENDIF
[3274]	612	IF ( bulk_cloud_model .AND. microphysics_seifert ) THEN
[1409]	613	qr_p(:,:,nxl-1) = qr_p(:,:,nxl)
	614	nr_p(:,:,nxl-1) = nr_p(:,:,nxl)
[1053]	615	ENDIF
[1409]	616	ENDIF
[1960]	617	IF ( passive_scalar ) s_p(:,:,nxl-1) = s_p(:,:,nxl)
[3182]	618	ELSEIF ( bc_radiation_r ) THEN
[1409]	619	pt_p(:,:,nxr+1) = pt_p(:,:,nxr)
[1960]	620	IF ( humidity ) THEN
[1409]	621	q_p(:,:,nxr+1) = q_p(:,:,nxr)
[3274]	622	IF ( bulk_cloud_model .AND. microphysics_morrison ) THEN
[2292]	623	qc_p(:,:,nxr+1) = qc_p(:,:,nxr)
	624	nc_p(:,:,nxr+1) = nc_p(:,:,nxr)
	625	ENDIF
[3274]	626	IF ( bulk_cloud_model .AND. microphysics_seifert ) THEN
[1409]	627	qr_p(:,:,nxr+1) = qr_p(:,:,nxr)
	628	nr_p(:,:,nxr+1) = nr_p(:,:,nxr)
[1053]	629	ENDIF
[1]	630	ENDIF
[1960]	631	IF ( passive_scalar ) s_p(:,:,nxr+1) = s_p(:,:,nxr)
[1]	632	ENDIF
	633
	634	!
[2696]	635	!-- Lateral boundary conditions for chemical species
	636	IF ( air_chemistry ) CALL chem_boundary_conds( 'set_bc_lateral' )
	637
	638	!
[1159]	639	!-- Radiation boundary conditions for the velocities at the respective outflow.
	640	!-- The phase velocity is either assumed to the maximum phase velocity that
	641	!-- ensures numerical stability (CFL-condition) or calculated after
	642	!-- Orlanski(1976) and averaged along the outflow boundary.
[3182]	643	IF ( bc_radiation_s ) THEN
[75]	644
[1159]	645	IF ( use_cmax ) THEN
	646	u_p(:,-1,:) = u(:,0,:)
	647	v_p(:,0,:) = v(:,1,:)
	648	w_p(:,-1,:) = w(:,0,:)
	649	ELSEIF ( .NOT. use_cmax ) THEN
[75]	650
[978]	651	c_max = dy / dt_3d
[75]	652
[1353]	653	c_u_m_l = 0.0_wp
	654	c_v_m_l = 0.0_wp
	655	c_w_m_l = 0.0_wp
[978]	656
[1353]	657	c_u_m = 0.0_wp
	658	c_v_m = 0.0_wp
	659	c_w_m = 0.0_wp
[978]	660
[75]	661	!
[996]	662	!-- Calculate the phase speeds for u, v, and w, first local and then
	663	!-- average along the outflow boundary.
	664	DO k = nzb+1, nzt+1
	665	DO i = nxl, nxr
[75]	666
[106]	667	denom = u_m_s(k,0,i) - u_m_s(k,1,i)
	668
[1353]	669	IF ( denom /= 0.0_wp ) THEN
[996]	670	c_u(k,i) = -c_max * ( u(k,0,i) - u_m_s(k,0,i) ) / ( denom * tsc(2) )
[1353]	671	IF ( c_u(k,i) < 0.0_wp ) THEN
	672	c_u(k,i) = 0.0_wp
[106]	673	ELSEIF ( c_u(k,i) > c_max ) THEN
	674	c_u(k,i) = c_max
	675	ENDIF
	676	ELSE
	677	c_u(k,i) = c_max
[75]	678	ENDIF
	679
[106]	680	denom = v_m_s(k,1,i) - v_m_s(k,2,i)
	681
[1353]	682	IF ( denom /= 0.0_wp ) THEN
[996]	683	c_v(k,i) = -c_max * ( v(k,1,i) - v_m_s(k,1,i) ) / ( denom * tsc(2) )
[1353]	684	IF ( c_v(k,i) < 0.0_wp ) THEN
	685	c_v(k,i) = 0.0_wp
[106]	686	ELSEIF ( c_v(k,i) > c_max ) THEN
	687	c_v(k,i) = c_max
	688	ENDIF
	689	ELSE
	690	c_v(k,i) = c_max
[75]	691	ENDIF
	692
[106]	693	denom = w_m_s(k,0,i) - w_m_s(k,1,i)
[75]	694
[1353]	695	IF ( denom /= 0.0_wp ) THEN
[996]	696	c_w(k,i) = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) / ( denom * tsc(2) )
[1353]	697	IF ( c_w(k,i) < 0.0_wp ) THEN
	698	c_w(k,i) = 0.0_wp
[106]	699	ELSEIF ( c_w(k,i) > c_max ) THEN
	700	c_w(k,i) = c_max
	701	ENDIF
	702	ELSE
	703	c_w(k,i) = c_max
[75]	704	ENDIF
[106]	705
[978]	706	c_u_m_l(k) = c_u_m_l(k) + c_u(k,i)
	707	c_v_m_l(k) = c_v_m_l(k) + c_v(k,i)
	708	c_w_m_l(k) = c_w_m_l(k) + c_w(k,i)
[106]	709
[978]	710	ENDDO
	711	ENDDO
[75]	712
[978]	713	#if defined( __parallel )
	714	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
	715	CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
	716	MPI_SUM, comm1dx, ierr )
	717	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
	718	CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
	719	MPI_SUM, comm1dx, ierr )
	720	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
	721	CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
	722	MPI_SUM, comm1dx, ierr )
	723	#else
	724	c_u_m = c_u_m_l
	725	c_v_m = c_v_m_l
	726	c_w_m = c_w_m_l
	727	#endif
	728
	729	c_u_m = c_u_m / (nx+1)
	730	c_v_m = c_v_m / (nx+1)
	731	c_w_m = c_w_m / (nx+1)
	732
[75]	733	!
[978]	734	!-- Save old timelevels for the next timestep
	735	IF ( intermediate_timestep_count == 1 ) THEN
	736	u_m_s(:,:,:) = u(:,0:1,:)
	737	v_m_s(:,:,:) = v(:,1:2,:)
	738	w_m_s(:,:,:) = w(:,0:1,:)
	739	ENDIF
	740
	741	!
	742	!-- Calculate the new velocities
[996]	743	DO k = nzb+1, nzt+1
	744	DO i = nxlg, nxrg
[978]	745	u_p(k,-1,i) = u(k,-1,i) - dt_3d * tsc(2) * c_u_m(k) * &
[75]	746	( u(k,-1,i) - u(k,0,i) ) * ddy
	747
[978]	748	v_p(k,0,i) = v(k,0,i) - dt_3d * tsc(2) * c_v_m(k) * &
[106]	749	( v(k,0,i) - v(k,1,i) ) * ddy
[75]	750
[978]	751	w_p(k,-1,i) = w(k,-1,i) - dt_3d * tsc(2) * c_w_m(k) * &
[75]	752	( w(k,-1,i) - w(k,0,i) ) * ddy
[978]	753	ENDDO
[75]	754	ENDDO
	755
	756	!
[978]	757	!-- Bottom boundary at the outflow
	758	IF ( ibc_uv_b == 0 ) THEN
[1353]	759	u_p(nzb,-1,:) = 0.0_wp
	760	v_p(nzb,0,:) = 0.0_wp
[978]	761	ELSE
	762	u_p(nzb,-1,:) = u_p(nzb+1,-1,:)
	763	v_p(nzb,0,:) = v_p(nzb+1,0,:)
	764	ENDIF
[1353]	765	w_p(nzb,-1,:) = 0.0_wp
[73]	766
[75]	767	!
[978]	768	!-- Top boundary at the outflow
	769	IF ( ibc_uv_t == 0 ) THEN
	770	u_p(nzt+1,-1,:) = u_init(nzt+1)
	771	v_p(nzt+1,0,:) = v_init(nzt+1)
	772	ELSE
[1742]	773	u_p(nzt+1,-1,:) = u_p(nzt,-1,:)
	774	v_p(nzt+1,0,:) = v_p(nzt,0,:)
[978]	775	ENDIF
[1353]	776	w_p(nzt:nzt+1,-1,:) = 0.0_wp
[978]	777
[75]	778	ENDIF
[73]	779
[75]	780	ENDIF
[73]	781
[3182]	782	IF ( bc_radiation_n ) THEN
[73]	783
[1159]	784	IF ( use_cmax ) THEN
	785	u_p(:,ny+1,:) = u(:,ny,:)
	786	v_p(:,ny+1,:) = v(:,ny,:)
	787	w_p(:,ny+1,:) = w(:,ny,:)
	788	ELSEIF ( .NOT. use_cmax ) THEN
[75]	789
[978]	790	c_max = dy / dt_3d
[75]	791
[1353]	792	c_u_m_l = 0.0_wp
	793	c_v_m_l = 0.0_wp
	794	c_w_m_l = 0.0_wp
[978]	795
[1353]	796	c_u_m = 0.0_wp
	797	c_v_m = 0.0_wp
	798	c_w_m = 0.0_wp
[978]	799
[1]	800	!
[996]	801	!-- Calculate the phase speeds for u, v, and w, first local and then
	802	!-- average along the outflow boundary.
	803	DO k = nzb+1, nzt+1
	804	DO i = nxl, nxr
[73]	805
[106]	806	denom = u_m_n(k,ny,i) - u_m_n(k,ny-1,i)
	807
[1353]	808	IF ( denom /= 0.0_wp ) THEN
[996]	809	c_u(k,i) = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) / ( denom * tsc(2) )
[1353]	810	IF ( c_u(k,i) < 0.0_wp ) THEN
	811	c_u(k,i) = 0.0_wp
[106]	812	ELSEIF ( c_u(k,i) > c_max ) THEN
	813	c_u(k,i) = c_max
	814	ENDIF
	815	ELSE
	816	c_u(k,i) = c_max
[73]	817	ENDIF
	818
[106]	819	denom = v_m_n(k,ny,i) - v_m_n(k,ny-1,i)
[73]	820
[1353]	821	IF ( denom /= 0.0_wp ) THEN
[996]	822	c_v(k,i) = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / ( denom * tsc(2) )
[1353]	823	IF ( c_v(k,i) < 0.0_wp ) THEN
	824	c_v(k,i) = 0.0_wp
[106]	825	ELSEIF ( c_v(k,i) > c_max ) THEN
	826	c_v(k,i) = c_max
	827	ENDIF
	828	ELSE
	829	c_v(k,i) = c_max
[73]	830	ENDIF
	831
[106]	832	denom = w_m_n(k,ny,i) - w_m_n(k,ny-1,i)
[73]	833
[1353]	834	IF ( denom /= 0.0_wp ) THEN
[996]	835	c_w(k,i) = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / ( denom * tsc(2) )
[1353]	836	IF ( c_w(k,i) < 0.0_wp ) THEN
	837	c_w(k,i) = 0.0_wp
[106]	838	ELSEIF ( c_w(k,i) > c_max ) THEN
	839	c_w(k,i) = c_max
	840	ENDIF
	841	ELSE
	842	c_w(k,i) = c_max
[73]	843	ENDIF
[106]	844
[978]	845	c_u_m_l(k) = c_u_m_l(k) + c_u(k,i)
	846	c_v_m_l(k) = c_v_m_l(k) + c_v(k,i)
	847	c_w_m_l(k) = c_w_m_l(k) + c_w(k,i)
[106]	848
[978]	849	ENDDO
	850	ENDDO
[73]	851
[978]	852	#if defined( __parallel )
	853	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
	854	CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
	855	MPI_SUM, comm1dx, ierr )
	856	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
	857	CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
	858	MPI_SUM, comm1dx, ierr )
	859	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
	860	CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
	861	MPI_SUM, comm1dx, ierr )
	862	#else
	863	c_u_m = c_u_m_l
	864	c_v_m = c_v_m_l
	865	c_w_m = c_w_m_l
	866	#endif
	867
	868	c_u_m = c_u_m / (nx+1)
	869	c_v_m = c_v_m / (nx+1)
	870	c_w_m = c_w_m / (nx+1)
	871
[73]	872	!
[978]	873	!-- Save old timelevels for the next timestep
	874	IF ( intermediate_timestep_count == 1 ) THEN
	875	u_m_n(:,:,:) = u(:,ny-1:ny,:)
	876	v_m_n(:,:,:) = v(:,ny-1:ny,:)
	877	w_m_n(:,:,:) = w(:,ny-1:ny,:)
	878	ENDIF
[73]	879
[978]	880	!
	881	!-- Calculate the new velocities
[996]	882	DO k = nzb+1, nzt+1
	883	DO i = nxlg, nxrg
[978]	884	u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * tsc(2) * c_u_m(k) * &
	885	( u(k,ny+1,i) - u(k,ny,i) ) * ddy
[73]	886
[978]	887	v_p(k,ny+1,i) = v(k,ny+1,i) - dt_3d * tsc(2) * c_v_m(k) * &
	888	( v(k,ny+1,i) - v(k,ny,i) ) * ddy
[73]	889
[978]	890	w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * tsc(2) * c_w_m(k) * &
	891	( w(k,ny+1,i) - w(k,ny,i) ) * ddy
	892	ENDDO
[1]	893	ENDDO
	894
	895	!
[978]	896	!-- Bottom boundary at the outflow
	897	IF ( ibc_uv_b == 0 ) THEN
[1353]	898	u_p(nzb,ny+1,:) = 0.0_wp
	899	v_p(nzb,ny+1,:) = 0.0_wp
[978]	900	ELSE
	901	u_p(nzb,ny+1,:) = u_p(nzb+1,ny+1,:)
	902	v_p(nzb,ny+1,:) = v_p(nzb+1,ny+1,:)
	903	ENDIF
[1353]	904	w_p(nzb,ny+1,:) = 0.0_wp
[73]	905
	906	!
[978]	907	!-- Top boundary at the outflow
	908	IF ( ibc_uv_t == 0 ) THEN
	909	u_p(nzt+1,ny+1,:) = u_init(nzt+1)
	910	v_p(nzt+1,ny+1,:) = v_init(nzt+1)
	911	ELSE
	912	u_p(nzt+1,ny+1,:) = u_p(nzt,nyn+1,:)
	913	v_p(nzt+1,ny+1,:) = v_p(nzt,nyn+1,:)
	914	ENDIF
[1353]	915	w_p(nzt:nzt+1,ny+1,:) = 0.0_wp
[978]	916
[1]	917	ENDIF
	918
[75]	919	ENDIF
	920
[3182]	921	IF ( bc_radiation_l ) THEN
[75]	922
[1159]	923	IF ( use_cmax ) THEN
[1717]	924	u_p(:,:,0) = u(:,:,1)
	925	v_p(:,:,-1) = v(:,:,0)
[1159]	926	w_p(:,:,-1) = w(:,:,0)
	927	ELSEIF ( .NOT. use_cmax ) THEN
[75]	928
[978]	929	c_max = dx / dt_3d
[75]	930
[1353]	931	c_u_m_l = 0.0_wp
	932	c_v_m_l = 0.0_wp
	933	c_w_m_l = 0.0_wp
[978]	934
[1353]	935	c_u_m = 0.0_wp
	936	c_v_m = 0.0_wp
	937	c_w_m = 0.0_wp
[978]	938
[1]	939	!
[996]	940	!-- Calculate the phase speeds for u, v, and w, first local and then
	941	!-- average along the outflow boundary.
	942	DO k = nzb+1, nzt+1
	943	DO j = nys, nyn
[75]	944
[106]	945	denom = u_m_l(k,j,1) - u_m_l(k,j,2)
	946
[1353]	947	IF ( denom /= 0.0_wp ) THEN
[996]	948	c_u(k,j) = -c_max * ( u(k,j,1) - u_m_l(k,j,1) ) / ( denom * tsc(2) )
[1353]	949	IF ( c_u(k,j) < 0.0_wp ) THEN
	950	c_u(k,j) = 0.0_wp
[107]	951	ELSEIF ( c_u(k,j) > c_max ) THEN
	952	c_u(k,j) = c_max
[106]	953	ENDIF
	954	ELSE
[107]	955	c_u(k,j) = c_max
[75]	956	ENDIF
	957
[106]	958	denom = v_m_l(k,j,0) - v_m_l(k,j,1)
[75]	959
[1353]	960	IF ( denom /= 0.0_wp ) THEN
[996]	961	c_v(k,j) = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) / ( denom * tsc(2) )
[1353]	962	IF ( c_v(k,j) < 0.0_wp ) THEN
	963	c_v(k,j) = 0.0_wp
[106]	964	ELSEIF ( c_v(k,j) > c_max ) THEN
	965	c_v(k,j) = c_max
	966	ENDIF
	967	ELSE
	968	c_v(k,j) = c_max
[75]	969	ENDIF
	970
[106]	971	denom = w_m_l(k,j,0) - w_m_l(k,j,1)
[75]	972
[1353]	973	IF ( denom /= 0.0_wp ) THEN
[996]	974	c_w(k,j) = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) / ( denom * tsc(2) )
[1353]	975	IF ( c_w(k,j) < 0.0_wp ) THEN
	976	c_w(k,j) = 0.0_wp
[106]	977	ELSEIF ( c_w(k,j) > c_max ) THEN
	978	c_w(k,j) = c_max
	979	ENDIF
	980	ELSE
	981	c_w(k,j) = c_max
[75]	982	ENDIF
[106]	983
[978]	984	c_u_m_l(k) = c_u_m_l(k) + c_u(k,j)
	985	c_v_m_l(k) = c_v_m_l(k) + c_v(k,j)
	986	c_w_m_l(k) = c_w_m_l(k) + c_w(k,j)
[106]	987
[978]	988	ENDDO
	989	ENDDO
[75]	990
[978]	991	#if defined( __parallel )
	992	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
	993	CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
	994	MPI_SUM, comm1dy, ierr )
	995	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
	996	CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
	997	MPI_SUM, comm1dy, ierr )
	998	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
	999	CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
	1000	MPI_SUM, comm1dy, ierr )
	1001	#else
	1002	c_u_m = c_u_m_l
	1003	c_v_m = c_v_m_l
	1004	c_w_m = c_w_m_l
	1005	#endif
	1006
	1007	c_u_m = c_u_m / (ny+1)
	1008	c_v_m = c_v_m / (ny+1)
	1009	c_w_m = c_w_m / (ny+1)
	1010
[73]	1011	!
[978]	1012	!-- Save old timelevels for the next timestep
	1013	IF ( intermediate_timestep_count == 1 ) THEN
	1014	u_m_l(:,:,:) = u(:,:,1:2)
	1015	v_m_l(:,:,:) = v(:,:,0:1)
	1016	w_m_l(:,:,:) = w(:,:,0:1)
	1017	ENDIF
	1018
	1019	!
	1020	!-- Calculate the new velocities
[996]	1021	DO k = nzb+1, nzt+1
[1113]	1022	DO j = nysg, nyng
[978]	1023	u_p(k,j,0) = u(k,j,0) - dt_3d * tsc(2) * c_u_m(k) * &
[106]	1024	( u(k,j,0) - u(k,j,1) ) * ddx
[75]	1025
[978]	1026	v_p(k,j,-1) = v(k,j,-1) - dt_3d * tsc(2) * c_v_m(k) * &
[75]	1027	( v(k,j,-1) - v(k,j,0) ) * ddx
	1028
[978]	1029	w_p(k,j,-1) = w(k,j,-1) - dt_3d * tsc(2) * c_w_m(k) * &
[75]	1030	( w(k,j,-1) - w(k,j,0) ) * ddx
[978]	1031	ENDDO
[75]	1032	ENDDO
	1033
	1034	!
[978]	1035	!-- Bottom boundary at the outflow
	1036	IF ( ibc_uv_b == 0 ) THEN
[1353]	1037	u_p(nzb,:,0) = 0.0_wp
	1038	v_p(nzb,:,-1) = 0.0_wp
[978]	1039	ELSE
	1040	u_p(nzb,:,0) = u_p(nzb+1,:,0)
	1041	v_p(nzb,:,-1) = v_p(nzb+1,:,-1)
	1042	ENDIF
[1353]	1043	w_p(nzb,:,-1) = 0.0_wp
[1]	1044
[75]	1045	!
[978]	1046	!-- Top boundary at the outflow
	1047	IF ( ibc_uv_t == 0 ) THEN
[1764]	1048	u_p(nzt+1,:,0) = u_init(nzt+1)
[978]	1049	v_p(nzt+1,:,-1) = v_init(nzt+1)
	1050	ELSE
[1764]	1051	u_p(nzt+1,:,0) = u_p(nzt,:,0)
[978]	1052	v_p(nzt+1,:,-1) = v_p(nzt,:,-1)
	1053	ENDIF
[1353]	1054	w_p(nzt:nzt+1,:,-1) = 0.0_wp
[978]	1055
[75]	1056	ENDIF
[73]	1057
[75]	1058	ENDIF
[73]	1059
[3182]	1060	IF ( bc_radiation_r ) THEN
[73]	1061
[1159]	1062	IF ( use_cmax ) THEN
	1063	u_p(:,:,nx+1) = u(:,:,nx)
	1064	v_p(:,:,nx+1) = v(:,:,nx)
	1065	w_p(:,:,nx+1) = w(:,:,nx)
	1066	ELSEIF ( .NOT. use_cmax ) THEN
[75]	1067
[978]	1068	c_max = dx / dt_3d
[75]	1069
[1353]	1070	c_u_m_l = 0.0_wp
	1071	c_v_m_l = 0.0_wp
	1072	c_w_m_l = 0.0_wp
[978]	1073
[1353]	1074	c_u_m = 0.0_wp
	1075	c_v_m = 0.0_wp
	1076	c_w_m = 0.0_wp
[978]	1077
[1]	1078	!
[996]	1079	!-- Calculate the phase speeds for u, v, and w, first local and then
	1080	!-- average along the outflow boundary.
	1081	DO k = nzb+1, nzt+1
	1082	DO j = nys, nyn
[73]	1083
[106]	1084	denom = u_m_r(k,j,nx) - u_m_r(k,j,nx-1)
	1085
[1353]	1086	IF ( denom /= 0.0_wp ) THEN
[996]	1087	c_u(k,j) = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) ) / ( denom * tsc(2) )
[1353]	1088	IF ( c_u(k,j) < 0.0_wp ) THEN
	1089	c_u(k,j) = 0.0_wp
[106]	1090	ELSEIF ( c_u(k,j) > c_max ) THEN
	1091	c_u(k,j) = c_max
	1092	ENDIF
	1093	ELSE
	1094	c_u(k,j) = c_max
[73]	1095	ENDIF
	1096
[106]	1097	denom = v_m_r(k,j,nx) - v_m_r(k,j,nx-1)
[73]	1098
[1353]	1099	IF ( denom /= 0.0_wp ) THEN
[996]	1100	c_v(k,j) = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) / ( denom * tsc(2) )
[1353]	1101	IF ( c_v(k,j) < 0.0_wp ) THEN
	1102	c_v(k,j) = 0.0_wp
[106]	1103	ELSEIF ( c_v(k,j) > c_max ) THEN
	1104	c_v(k,j) = c_max
	1105	ENDIF
	1106	ELSE
	1107	c_v(k,j) = c_max
[73]	1108	ENDIF
	1109
[106]	1110	denom = w_m_r(k,j,nx) - w_m_r(k,j,nx-1)
[73]	1111
[1353]	1112	IF ( denom /= 0.0_wp ) THEN
[996]	1113	c_w(k,j) = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) / ( denom * tsc(2) )
[1353]	1114	IF ( c_w(k,j) < 0.0_wp ) THEN
	1115	c_w(k,j) = 0.0_wp
[106]	1116	ELSEIF ( c_w(k,j) > c_max ) THEN
	1117	c_w(k,j) = c_max
	1118	ENDIF
	1119	ELSE
	1120	c_w(k,j) = c_max
[73]	1121	ENDIF
[106]	1122
[978]	1123	c_u_m_l(k) = c_u_m_l(k) + c_u(k,j)
	1124	c_v_m_l(k) = c_v_m_l(k) + c_v(k,j)
	1125	c_w_m_l(k) = c_w_m_l(k) + c_w(k,j)
[106]	1126
[978]	1127	ENDDO
	1128	ENDDO
[73]	1129
[978]	1130	#if defined( __parallel )
	1131	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
	1132	CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
	1133	MPI_SUM, comm1dy, ierr )
	1134	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
	1135	CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
	1136	MPI_SUM, comm1dy, ierr )
	1137	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
	1138	CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
	1139	MPI_SUM, comm1dy, ierr )
	1140	#else
	1141	c_u_m = c_u_m_l
	1142	c_v_m = c_v_m_l
	1143	c_w_m = c_w_m_l
	1144	#endif
	1145
	1146	c_u_m = c_u_m / (ny+1)
	1147	c_v_m = c_v_m / (ny+1)
	1148	c_w_m = c_w_m / (ny+1)
	1149
[73]	1150	!
[978]	1151	!-- Save old timelevels for the next timestep
	1152	IF ( intermediate_timestep_count == 1 ) THEN
	1153	u_m_r(:,:,:) = u(:,:,nx-1:nx)
	1154	v_m_r(:,:,:) = v(:,:,nx-1:nx)
	1155	w_m_r(:,:,:) = w(:,:,nx-1:nx)
	1156	ENDIF
[73]	1157
[978]	1158	!
	1159	!-- Calculate the new velocities
[996]	1160	DO k = nzb+1, nzt+1
[1113]	1161	DO j = nysg, nyng
[978]	1162	u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * tsc(2) * c_u_m(k) * &
	1163	( u(k,j,nx+1) - u(k,j,nx) ) * ddx
[73]	1164
[978]	1165	v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * tsc(2) * c_v_m(k) * &
	1166	( v(k,j,nx+1) - v(k,j,nx) ) * ddx
[73]	1167
[978]	1168	w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * tsc(2) * c_w_m(k) * &
	1169	( w(k,j,nx+1) - w(k,j,nx) ) * ddx
	1170	ENDDO
[73]	1171	ENDDO
	1172
	1173	!
[978]	1174	!-- Bottom boundary at the outflow
	1175	IF ( ibc_uv_b == 0 ) THEN
[1353]	1176	u_p(nzb,:,nx+1) = 0.0_wp
	1177	v_p(nzb,:,nx+1) = 0.0_wp
[978]	1178	ELSE
	1179	u_p(nzb,:,nx+1) = u_p(nzb+1,:,nx+1)
	1180	v_p(nzb,:,nx+1) = v_p(nzb+1,:,nx+1)
	1181	ENDIF
[1353]	1182	w_p(nzb,:,nx+1) = 0.0_wp
[73]	1183
	1184	!
[978]	1185	!-- Top boundary at the outflow
	1186	IF ( ibc_uv_t == 0 ) THEN
	1187	u_p(nzt+1,:,nx+1) = u_init(nzt+1)
	1188	v_p(nzt+1,:,nx+1) = v_init(nzt+1)
	1189	ELSE
	1190	u_p(nzt+1,:,nx+1) = u_p(nzt,:,nx+1)
	1191	v_p(nzt+1,:,nx+1) = v_p(nzt,:,nx+1)
	1192	ENDIF
[1742]	1193	w_p(nzt:nzt+1,:,nx+1) = 0.0_wp
[978]	1194
[1]	1195	ENDIF
	1196
	1197	ENDIF
[3864]	1198
[3467]	1199	IF ( salsa ) THEN
	1200	CALL salsa_boundary_conds
[3864]	1201	ENDIF
[1]	1202
	1203	END SUBROUTINE boundary_conds

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |