Home

Context Navigation

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

Last change on this file since 4212 was 4182, checked in by scharf, 5 years ago
corrected "Former revisions" section minor formatting in "Former revisions" section added "Author" section
Property svn:keywords set to `Id`
File size: 35.0 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 4182 2019-08-22 15:20:23Z suehring $
[4182]	27	! Corrected "Former revisions" section
	28	!
	29	! 4102 2019-07-17 16:00:03Z suehring
[4102]	30	! - Revise setting for boundary conditions at horizontal walls, use the offset
	31	! index that belongs to the data structure instead of pre-calculating
	32	! the offset index for each facing.
	33	! - Set boundary conditions for bulk-cloud quantities also at downward facing
	34	! surfaces
	35	!
	36	! 4087 2019-07-11 11:35:23Z gronemeier
[4087]	37	! Add comment
	38	!
	39	! 4086 2019-07-11 05:55:44Z gronemeier
[4086]	40	! Bugfix: use constant-flux layer condition for e in all rans modes
	41	!
	42	! 3879 2019-04-08 20:25:23Z knoop
[3717]	43	! Bugfix, do not set boundary conditions for potential temperature in neutral
	44	! runs.
	45	!
	46	! 3655 2019-01-07 16:51:22Z knoop
[3634]	47	! OpenACC port for SPEC
[1321]	48	!
[4182]	49	! Revision 1.1 1997/09/12 06:21:34 raasch
	50	! Initial revision
	51	!
	52	!
[1]	53	! Description:
	54	! ------------
[1682]	55	!> Boundary conditions for the prognostic quantities.
	56	!> One additional bottom boundary condition is applied for the TKE (=(u)*2)
	57	!> in prandtl_fluxes. The cyclic lateral boundary conditions are implicitly
	58	!> handled in routine exchange_horiz. Pressure boundary conditions are
	59	!> explicitly set in routines pres, poisfft, poismg and sor.
[1]	60	!------------------------------------------------------------------------------!
[1682]	61	SUBROUTINE boundary_conds
	62
[1]	63
[1320]	64	USE arrays_3d, &
	65	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]	66	dzu, nc_p, nr_p, pt, pt_init, pt_p, q, &
[3241]	67	q_p, qc_p, qr_p, s, s_p, sa, sa_p, u, u_init, u_m_l, u_m_n, &
	68	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, &
	69	w, w_p, w_m_l, w_m_n, w_m_r, w_m_s
[2696]	70
[3294]	71	USE bulk_cloud_model_mod, &
	72	ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert
	73
[2696]	74	USE chemistry_model_mod, &
[3879]	75	ONLY: chem_boundary_conds
	76
[1320]	77	USE control_parameters, &
[4102]	78	ONLY: air_chemistry, bc_dirichlet_l, &
[3182]	79	bc_dirichlet_s, bc_radiation_l, bc_radiation_n, bc_radiation_r, &
	80	bc_radiation_s, bc_pt_t_val, bc_q_t_val, bc_s_t_val, &
[4102]	81	child_domain, coupling_mode, dt_3d, &
[3582]	82	humidity, ibc_pt_b, ibc_pt_t, ibc_q_b, ibc_q_t, ibc_s_b, &
	83	ibc_s_t, ibc_uv_b, ibc_uv_t, intermediate_timestep_count, &
[3717]	84	nesting_offline, neutral, nudging, ocean_mode, passive_scalar, &
[4102]	85	tsc, salsa, use_cmax
[1320]	86
	87	USE grid_variables, &
	88	ONLY: ddx, ddy, dx, dy
	89
	90	USE indices, &
[3294]	91	ONLY: nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb, nzt
[1320]	92
	93	USE kinds
	94
[3294]	95	USE ocean_mod, &
	96	ONLY: ibc_sa_t
[3274]	97
[1]	98	USE pegrid
	99
[1933]	100	USE pmc_interface, &
[4102]	101	ONLY : nesting_mode
[3467]	102
	103	USE salsa_mod, &
[3582]	104	ONLY: salsa_boundary_conds
[1320]	105
[2232]	106	USE surface_mod, &
[4102]	107	ONLY : bc_h
[1933]	108
[3129]	109	USE turbulence_closure_mod, &
[4102]	110	ONLY: tcm_boundary_conds
[3129]	111
[1]	112	IMPLICIT NONE
	113
[2232]	114	INTEGER(iwp) :: i !< grid index x direction
	115	INTEGER(iwp) :: j !< grid index y direction
	116	INTEGER(iwp) :: k !< grid index z direction
	117	INTEGER(iwp) :: l !< running index boundary type, for up- and downward-facing walls
	118	INTEGER(iwp) :: m !< running index surface elements
[1]	119
[3562]	120	REAL(wp) :: c_max !< maximum phase velocity allowed by CFL criterion, used for outflow boundary condition
	121	REAL(wp) :: denom !< horizontal gradient of velocity component normal to the outflow boundary
[1]	122
[73]	123
[1]	124	!
[1113]	125	!-- Bottom boundary
	126	IF ( ibc_uv_b == 1 ) THEN
	127	u_p(nzb,:,:) = u_p(nzb+1,:,:)
	128	v_p(nzb,:,:) = v_p(nzb+1,:,:)
	129	ENDIF
[2232]	130	!
	131	!-- Set zero vertical velocity at topography top (l=0), or bottom (l=1) in case
	132	!-- of downward-facing surfaces.
	133	DO l = 0, 1
	134	!$OMP PARALLEL DO PRIVATE( i, j, k )
[3634]	135	!$ACC PARALLEL LOOP PRIVATE(i, j, k) &
	136	!$ACC PRESENT(bc_h, w_p)
[2232]	137	DO m = 1, bc_h(l)%ns
	138	i = bc_h(l)%i(m)
	139	j = bc_h(l)%j(m)
	140	k = bc_h(l)%k(m)
[4102]	141	w_p(k+bc_h(l)%koff,j,i) = 0.0_wp
[1113]	142	ENDDO
	143	ENDDO
	144
	145	!
[1762]	146	!-- Top boundary. A nested domain ( ibc_uv_t = 3 ) does not require settings.
[1113]	147	IF ( ibc_uv_t == 0 ) THEN
[3634]	148	!$ACC KERNELS PRESENT(u_p, v_p, u_init, v_init)
[1113]	149	u_p(nzt+1,:,:) = u_init(nzt+1)
	150	v_p(nzt+1,:,:) = v_init(nzt+1)
[3634]	151	!$ACC END KERNELS
[1762]	152	ELSEIF ( ibc_uv_t == 1 ) THEN
[1113]	153	u_p(nzt+1,:,:) = u_p(nzt,:,:)
	154	v_p(nzt+1,:,:) = v_p(nzt,:,:)
	155	ENDIF
	156
[2365]	157	!
	158	!-- Vertical nesting: Vertical velocity not zero at the top of the fine grid
[3347]	159	IF ( .NOT. child_domain .AND. .NOT. nesting_offline .AND. &
[2365]	160	TRIM(coupling_mode) /= 'vnested_fine' ) THEN
[3634]	161	!$ACC KERNELS PRESENT(w_p)
[2365]	162	w_p(nzt:nzt+1,:,:) = 0.0_wp !< nzt is not a prognostic level (but cf. pres)
[3634]	163	!$ACC END KERNELS
[1762]	164	ENDIF
	165
[1113]	166	!
[2232]	167	!-- Temperature at bottom and top boundary.
[1113]	168	!-- In case of coupled runs (ibc_pt_b = 2) the temperature is given by
	169	!-- the sea surface temperature of the coupled ocean model.
[2232]	170	!-- Dirichlet
[3717]	171	IF ( .NOT. neutral ) THEN
	172	IF ( ibc_pt_b == 0 ) THEN
	173	DO l = 0, 1
	174	!$OMP PARALLEL DO PRIVATE( i, j, k )
	175	DO m = 1, bc_h(l)%ns
	176	i = bc_h(l)%i(m)
	177	j = bc_h(l)%j(m)
	178	k = bc_h(l)%k(m)
[4102]	179	pt_p(k+bc_h(l)%koff,j,i) = pt(k+bc_h(l)%koff,j,i)
[3717]	180	ENDDO
[1]	181	ENDDO
[3717]	182	!
	183	!-- Neumann, zero-gradient
	184	ELSEIF ( ibc_pt_b == 1 ) THEN
	185	DO l = 0, 1
	186	!$OMP PARALLEL DO PRIVATE( i, j, k )
	187	!$ACC PARALLEL LOOP PRIVATE(i, j, k) &
	188	!$ACC PRESENT(bc_h, pt_p)
	189	DO m = 1, bc_h(l)%ns
	190	i = bc_h(l)%i(m)
	191	j = bc_h(l)%j(m)
	192	k = bc_h(l)%k(m)
[4102]	193	pt_p(k+bc_h(l)%koff,j,i) = pt_p(k,j,i)
[3717]	194	ENDDO
[1113]	195	ENDDO
[3717]	196	ENDIF
	197
	198	!
	199	!-- Temperature at top boundary
	200	IF ( ibc_pt_t == 0 ) THEN
	201	pt_p(nzt+1,:,:) = pt(nzt+1,:,:)
	202	!
	203	!-- In case of nudging adjust top boundary to pt which is
	204	!-- read in from NUDGING-DATA
	205	IF ( nudging ) THEN
	206	pt_p(nzt+1,:,:) = pt_init(nzt+1)
	207	ENDIF
	208	ELSEIF ( ibc_pt_t == 1 ) THEN
	209	pt_p(nzt+1,:,:) = pt_p(nzt,:,:)
	210	ELSEIF ( ibc_pt_t == 2 ) THEN
	211	!$ACC KERNELS PRESENT(pt_p, dzu)
	212	pt_p(nzt+1,:,:) = pt_p(nzt,:,:) + bc_pt_t_val * dzu(nzt+1)
	213	!$ACC END KERNELS
	214	ENDIF
[1113]	215	ENDIF
[1]	216	!
[1113]	217	!-- Boundary conditions for salinity
[3294]	218	IF ( ocean_mode ) THEN
[1113]	219	!
	220	!-- Bottom boundary: Neumann condition because salinity flux is always
[2232]	221	!-- given.
	222	DO l = 0, 1
	223	!$OMP PARALLEL DO PRIVATE( i, j, k )
	224	DO m = 1, bc_h(l)%ns
	225	i = bc_h(l)%i(m)
	226	j = bc_h(l)%j(m)
	227	k = bc_h(l)%k(m)
[4102]	228	sa_p(k+bc_h(l)%koff,j,i) = sa_p(k,j,i)
[1]	229	ENDDO
[1113]	230	ENDDO
[1]	231	!
[1113]	232	!-- Top boundary: Dirichlet or Neumann
	233	IF ( ibc_sa_t == 0 ) THEN
	234	sa_p(nzt+1,:,:) = sa(nzt+1,:,:)
	235	ELSEIF ( ibc_sa_t == 1 ) THEN
	236	sa_p(nzt+1,:,:) = sa_p(nzt,:,:)
[1]	237	ENDIF
	238
[1113]	239	ENDIF
	240
[1]	241	!
[1960]	242	!-- Boundary conditions for total water content,
[1113]	243	!-- bottom and top boundary (see also temperature)
[1960]	244	IF ( humidity ) THEN
[1113]	245	!
	246	!-- Surface conditions for constant_humidity_flux
[2232]	247	!-- Run loop over all non-natural and natural walls. Note, in wall-datatype
	248	!-- the k coordinate belongs to the atmospheric grid point, therefore, set
	249	!-- q_p at k-1
[1113]	250	IF ( ibc_q_b == 0 ) THEN
[2232]	251
	252	DO l = 0, 1
	253	!$OMP PARALLEL DO PRIVATE( i, j, k )
	254	DO m = 1, bc_h(l)%ns
	255	i = bc_h(l)%i(m)
	256	j = bc_h(l)%j(m)
	257	k = bc_h(l)%k(m)
[4102]	258	q_p(k+bc_h(l)%koff,j,i) = q(k+bc_h(l)%koff,j,i)
[1]	259	ENDDO
	260	ENDDO
[2232]	261
[1113]	262	ELSE
[2232]	263
	264	DO l = 0, 1
	265	!$OMP PARALLEL DO PRIVATE( i, j, k )
	266	DO m = 1, bc_h(l)%ns
	267	i = bc_h(l)%i(m)
	268	j = bc_h(l)%j(m)
	269	k = bc_h(l)%k(m)
[4102]	270	q_p(k+bc_h(l)%koff,j,i) = q_p(k,j,i)
[95]	271	ENDDO
	272	ENDDO
[1113]	273	ENDIF
[95]	274	!
[1113]	275	!-- Top boundary
[1462]	276	IF ( ibc_q_t == 0 ) THEN
	277	q_p(nzt+1,:,:) = q(nzt+1,:,:)
	278	ELSEIF ( ibc_q_t == 1 ) THEN
[1992]	279	q_p(nzt+1,:,:) = q_p(nzt,:,:) + bc_q_t_val * dzu(nzt+1)
[1462]	280	ENDIF
[95]	281
[3274]	282	IF ( bulk_cloud_model .AND. microphysics_morrison ) THEN
[2292]	283	!
	284	!-- Surface conditions cloud water (Dirichlet)
	285	!-- Run loop over all non-natural and natural walls. Note, in wall-datatype
	286	!-- the k coordinate belongs to the atmospheric grid point, therefore, set
[4102]	287	!-- qr_p and nr_p at upward (k-1) and downward-facing (k+1) walls
	288	DO l = 0, 1
[2292]	289	!$OMP PARALLEL DO PRIVATE( i, j, k )
[4102]	290	DO m = 1, bc_h(l)%ns
	291	i = bc_h(l)%i(m)
	292	j = bc_h(l)%j(m)
	293	k = bc_h(l)%k(m)
	294	qc_p(k+bc_h(l)%koff,j,i) = 0.0_wp
	295	nc_p(k+bc_h(l)%koff,j,i) = 0.0_wp
	296	ENDDO
[2292]	297	ENDDO
	298	!
	299	!-- Top boundary condition for cloud water (Dirichlet)
	300	qc_p(nzt+1,:,:) = 0.0_wp
	301	nc_p(nzt+1,:,:) = 0.0_wp
	302
	303	ENDIF
	304
[3274]	305	IF ( bulk_cloud_model .AND. microphysics_seifert ) THEN
[1113]	306	!
[1361]	307	!-- Surface conditions rain water (Dirichlet)
[2232]	308	!-- Run loop over all non-natural and natural walls. Note, in wall-datatype
	309	!-- the k coordinate belongs to the atmospheric grid point, therefore, set
[4102]	310	!-- qr_p and nr_p at upward (k-1) and downward-facing (k+1) walls
	311	DO l = 0, 1
[2232]	312	!$OMP PARALLEL DO PRIVATE( i, j, k )
[4102]	313	DO m = 1, bc_h(l)%ns
	314	i = bc_h(l)%i(m)
	315	j = bc_h(l)%j(m)
	316	k = bc_h(l)%k(m)
	317	qr_p(k+bc_h(l)%koff,j,i) = 0.0_wp
	318	nr_p(k+bc_h(l)%koff,j,i) = 0.0_wp
	319	ENDDO
[1115]	320	ENDDO
[1]	321	!
[1361]	322	!-- Top boundary condition for rain water (Dirichlet)
	323	qr_p(nzt+1,:,:) = 0.0_wp
	324	nr_p(nzt+1,:,:) = 0.0_wp
[1115]	325
[1]	326	ENDIF
[1409]	327	ENDIF
[1]	328	!
[1960]	329	!-- Boundary conditions for scalar,
	330	!-- bottom and top boundary (see also temperature)
	331	IF ( passive_scalar ) THEN
	332	!
	333	!-- Surface conditions for constant_humidity_flux
[2232]	334	!-- Run loop over all non-natural and natural walls. Note, in wall-datatype
	335	!-- the k coordinate belongs to the atmospheric grid point, therefore, set
	336	!-- s_p at k-1
[1960]	337	IF ( ibc_s_b == 0 ) THEN
[2232]	338
	339	DO l = 0, 1
	340	!$OMP PARALLEL DO PRIVATE( i, j, k )
	341	DO m = 1, bc_h(l)%ns
	342	i = bc_h(l)%i(m)
	343	j = bc_h(l)%j(m)
	344	k = bc_h(l)%k(m)
[4102]	345	s_p(k+bc_h(l)%koff,j,i) = s(k+bc_h(l)%koff,j,i)
[1960]	346	ENDDO
	347	ENDDO
[2232]	348
[1960]	349	ELSE
[2232]	350
	351	DO l = 0, 1
	352	!$OMP PARALLEL DO PRIVATE( i, j, k )
	353	DO m = 1, bc_h(l)%ns
	354	i = bc_h(l)%i(m)
	355	j = bc_h(l)%j(m)
	356	k = bc_h(l)%k(m)
[4102]	357	s_p(k+bc_h(l)%koff,j,i) = s_p(k,j,i)
[1960]	358	ENDDO
	359	ENDDO
	360	ENDIF
	361	!
[1992]	362	!-- Top boundary condition
	363	IF ( ibc_s_t == 0 ) THEN
[1960]	364	s_p(nzt+1,:,:) = s(nzt+1,:,:)
[1992]	365	ELSEIF ( ibc_s_t == 1 ) THEN
	366	s_p(nzt+1,:,:) = s_p(nzt,:,:)
	367	ELSEIF ( ibc_s_t == 2 ) THEN
	368	s_p(nzt+1,:,:) = s_p(nzt,:,:) + bc_s_t_val * dzu(nzt+1)
[1960]	369	ENDIF
	370
[4102]	371	ENDIF
[1960]	372	!
[4102]	373	!-- Set boundary conditions for subgrid TKE and dissipation (RANS mode)
	374	CALL tcm_boundary_conds
	375	!
[2696]	376	!-- Top/bottom boundary conditions for chemical species
	377	IF ( air_chemistry ) CALL chem_boundary_conds( 'set_bc_bottomtop' )
	378	!
[1762]	379	!-- In case of inflow or nest boundary at the south boundary the boundary for v
	380	!-- is at nys and in case of inflow or nest boundary at the left boundary the
	381	!-- boundary for u is at nxl. Since in prognostic_equations (cache optimized
	382	!-- version) these levels are handled as a prognostic level, boundary values
	383	!-- have to be restored here.
[3182]	384	IF ( bc_dirichlet_s ) THEN
[1409]	385	v_p(:,nys,:) = v_p(:,nys-1,:)
[3182]	386	ELSEIF ( bc_dirichlet_l ) THEN
[1409]	387	u_p(:,:,nxl) = u_p(:,:,nxl-1)
	388	ENDIF
[1]	389
	390	!
[1762]	391	!-- The same restoration for u at i=nxl and v at j=nys as above must be made
[1933]	392	!-- in case of nest boundaries. This must not be done in case of vertical nesting
[3182]	393	!-- mode as in that case the lateral boundaries are actually cyclic.
[4102]	394	!-- Lateral oundary conditions for TKE and dissipation are set
	395	!-- in tcm_boundary_conds.
[3182]	396	IF ( nesting_mode /= 'vertical' .OR. nesting_offline ) THEN
	397	IF ( bc_dirichlet_s ) THEN
[1933]	398	v_p(:,nys,:) = v_p(:,nys-1,:)
	399	ENDIF
[3182]	400	IF ( bc_dirichlet_l ) THEN
[1933]	401	u_p(:,:,nxl) = u_p(:,:,nxl-1)
	402	ENDIF
[1762]	403	ENDIF
	404
	405	!
[4102]	406	!-- Lateral boundary conditions for scalar quantities at the outflow.
	407	!-- Lateral oundary conditions for TKE and dissipation are set
	408	!-- in tcm_boundary_conds.
[3182]	409	IF ( bc_radiation_s ) THEN
[1409]	410	pt_p(:,nys-1,:) = pt_p(:,nys,:)
[1960]	411	IF ( humidity ) THEN
[1409]	412	q_p(:,nys-1,:) = q_p(:,nys,:)
[3274]	413	IF ( bulk_cloud_model .AND. microphysics_morrison ) THEN
[2292]	414	qc_p(:,nys-1,:) = qc_p(:,nys,:)
	415	nc_p(:,nys-1,:) = nc_p(:,nys,:)
	416	ENDIF
[3274]	417	IF ( bulk_cloud_model .AND. microphysics_seifert ) THEN
[1409]	418	qr_p(:,nys-1,:) = qr_p(:,nys,:)
	419	nr_p(:,nys-1,:) = nr_p(:,nys,:)
[1053]	420	ENDIF
[1409]	421	ENDIF
[1960]	422	IF ( passive_scalar ) s_p(:,nys-1,:) = s_p(:,nys,:)
[3182]	423	ELSEIF ( bc_radiation_n ) THEN
[1409]	424	pt_p(:,nyn+1,:) = pt_p(:,nyn,:)
[1960]	425	IF ( humidity ) THEN
[1409]	426	q_p(:,nyn+1,:) = q_p(:,nyn,:)
[3274]	427	IF ( bulk_cloud_model .AND. microphysics_morrison ) THEN
[2292]	428	qc_p(:,nyn+1,:) = qc_p(:,nyn,:)
	429	nc_p(:,nyn+1,:) = nc_p(:,nyn,:)
	430	ENDIF
[3274]	431	IF ( bulk_cloud_model .AND. microphysics_seifert ) THEN
[1409]	432	qr_p(:,nyn+1,:) = qr_p(:,nyn,:)
	433	nr_p(:,nyn+1,:) = nr_p(:,nyn,:)
[1053]	434	ENDIF
[1409]	435	ENDIF
[1960]	436	IF ( passive_scalar ) s_p(:,nyn+1,:) = s_p(:,nyn,:)
[3182]	437	ELSEIF ( bc_radiation_l ) THEN
[1409]	438	pt_p(:,:,nxl-1) = pt_p(:,:,nxl)
[1960]	439	IF ( humidity ) THEN
[1409]	440	q_p(:,:,nxl-1) = q_p(:,:,nxl)
[3274]	441	IF ( bulk_cloud_model .AND. microphysics_morrison ) THEN
[2292]	442	qc_p(:,:,nxl-1) = qc_p(:,:,nxl)
	443	nc_p(:,:,nxl-1) = nc_p(:,:,nxl)
	444	ENDIF
[3274]	445	IF ( bulk_cloud_model .AND. microphysics_seifert ) THEN
[1409]	446	qr_p(:,:,nxl-1) = qr_p(:,:,nxl)
	447	nr_p(:,:,nxl-1) = nr_p(:,:,nxl)
[1053]	448	ENDIF
[1409]	449	ENDIF
[1960]	450	IF ( passive_scalar ) s_p(:,:,nxl-1) = s_p(:,:,nxl)
[3182]	451	ELSEIF ( bc_radiation_r ) THEN
[1409]	452	pt_p(:,:,nxr+1) = pt_p(:,:,nxr)
[1960]	453	IF ( humidity ) THEN
[1409]	454	q_p(:,:,nxr+1) = q_p(:,:,nxr)
[3274]	455	IF ( bulk_cloud_model .AND. microphysics_morrison ) THEN
[2292]	456	qc_p(:,:,nxr+1) = qc_p(:,:,nxr)
	457	nc_p(:,:,nxr+1) = nc_p(:,:,nxr)
	458	ENDIF
[3274]	459	IF ( bulk_cloud_model .AND. microphysics_seifert ) THEN
[1409]	460	qr_p(:,:,nxr+1) = qr_p(:,:,nxr)
	461	nr_p(:,:,nxr+1) = nr_p(:,:,nxr)
[1053]	462	ENDIF
[1]	463	ENDIF
[1960]	464	IF ( passive_scalar ) s_p(:,:,nxr+1) = s_p(:,:,nxr)
[1]	465	ENDIF
	466
	467	!
[2696]	468	!-- Lateral boundary conditions for chemical species
	469	IF ( air_chemistry ) CALL chem_boundary_conds( 'set_bc_lateral' )
	470
	471	!
[1159]	472	!-- Radiation boundary conditions for the velocities at the respective outflow.
	473	!-- The phase velocity is either assumed to the maximum phase velocity that
	474	!-- ensures numerical stability (CFL-condition) or calculated after
	475	!-- Orlanski(1976) and averaged along the outflow boundary.
[3182]	476	IF ( bc_radiation_s ) THEN
[75]	477
[1159]	478	IF ( use_cmax ) THEN
	479	u_p(:,-1,:) = u(:,0,:)
	480	v_p(:,0,:) = v(:,1,:)
	481	w_p(:,-1,:) = w(:,0,:)
	482	ELSEIF ( .NOT. use_cmax ) THEN
[75]	483
[978]	484	c_max = dy / dt_3d
[75]	485
[1353]	486	c_u_m_l = 0.0_wp
	487	c_v_m_l = 0.0_wp
	488	c_w_m_l = 0.0_wp
[978]	489
[1353]	490	c_u_m = 0.0_wp
	491	c_v_m = 0.0_wp
	492	c_w_m = 0.0_wp
[978]	493
[75]	494	!
[996]	495	!-- Calculate the phase speeds for u, v, and w, first local and then
	496	!-- average along the outflow boundary.
	497	DO k = nzb+1, nzt+1
	498	DO i = nxl, nxr
[75]	499
[106]	500	denom = u_m_s(k,0,i) - u_m_s(k,1,i)
	501
[1353]	502	IF ( denom /= 0.0_wp ) THEN
[996]	503	c_u(k,i) = -c_max * ( u(k,0,i) - u_m_s(k,0,i) ) / ( denom * tsc(2) )
[1353]	504	IF ( c_u(k,i) < 0.0_wp ) THEN
	505	c_u(k,i) = 0.0_wp
[106]	506	ELSEIF ( c_u(k,i) > c_max ) THEN
	507	c_u(k,i) = c_max
	508	ENDIF
	509	ELSE
	510	c_u(k,i) = c_max
[75]	511	ENDIF
	512
[106]	513	denom = v_m_s(k,1,i) - v_m_s(k,2,i)
	514
[1353]	515	IF ( denom /= 0.0_wp ) THEN
[996]	516	c_v(k,i) = -c_max * ( v(k,1,i) - v_m_s(k,1,i) ) / ( denom * tsc(2) )
[1353]	517	IF ( c_v(k,i) < 0.0_wp ) THEN
	518	c_v(k,i) = 0.0_wp
[106]	519	ELSEIF ( c_v(k,i) > c_max ) THEN
	520	c_v(k,i) = c_max
	521	ENDIF
	522	ELSE
	523	c_v(k,i) = c_max
[75]	524	ENDIF
	525
[106]	526	denom = w_m_s(k,0,i) - w_m_s(k,1,i)
[75]	527
[1353]	528	IF ( denom /= 0.0_wp ) THEN
[996]	529	c_w(k,i) = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) / ( denom * tsc(2) )
[1353]	530	IF ( c_w(k,i) < 0.0_wp ) THEN
	531	c_w(k,i) = 0.0_wp
[106]	532	ELSEIF ( c_w(k,i) > c_max ) THEN
	533	c_w(k,i) = c_max
	534	ENDIF
	535	ELSE
	536	c_w(k,i) = c_max
[75]	537	ENDIF
[106]	538
[978]	539	c_u_m_l(k) = c_u_m_l(k) + c_u(k,i)
	540	c_v_m_l(k) = c_v_m_l(k) + c_v(k,i)
	541	c_w_m_l(k) = c_w_m_l(k) + c_w(k,i)
[106]	542
[978]	543	ENDDO
	544	ENDDO
[75]	545
[978]	546	#if defined( __parallel )
	547	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
	548	CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
	549	MPI_SUM, comm1dx, ierr )
	550	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
	551	CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
	552	MPI_SUM, comm1dx, ierr )
	553	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
	554	CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
	555	MPI_SUM, comm1dx, ierr )
	556	#else
	557	c_u_m = c_u_m_l
	558	c_v_m = c_v_m_l
	559	c_w_m = c_w_m_l
	560	#endif
	561
	562	c_u_m = c_u_m / (nx+1)
	563	c_v_m = c_v_m / (nx+1)
	564	c_w_m = c_w_m / (nx+1)
	565
[75]	566	!
[978]	567	!-- Save old timelevels for the next timestep
	568	IF ( intermediate_timestep_count == 1 ) THEN
	569	u_m_s(:,:,:) = u(:,0:1,:)
	570	v_m_s(:,:,:) = v(:,1:2,:)
	571	w_m_s(:,:,:) = w(:,0:1,:)
	572	ENDIF
	573
	574	!
	575	!-- Calculate the new velocities
[996]	576	DO k = nzb+1, nzt+1
	577	DO i = nxlg, nxrg
[978]	578	u_p(k,-1,i) = u(k,-1,i) - dt_3d * tsc(2) * c_u_m(k) * &
[75]	579	( u(k,-1,i) - u(k,0,i) ) * ddy
	580
[978]	581	v_p(k,0,i) = v(k,0,i) - dt_3d * tsc(2) * c_v_m(k) * &
[106]	582	( v(k,0,i) - v(k,1,i) ) * ddy
[75]	583
[978]	584	w_p(k,-1,i) = w(k,-1,i) - dt_3d * tsc(2) * c_w_m(k) * &
[75]	585	( w(k,-1,i) - w(k,0,i) ) * ddy
[978]	586	ENDDO
[75]	587	ENDDO
	588
	589	!
[978]	590	!-- Bottom boundary at the outflow
	591	IF ( ibc_uv_b == 0 ) THEN
[1353]	592	u_p(nzb,-1,:) = 0.0_wp
	593	v_p(nzb,0,:) = 0.0_wp
[978]	594	ELSE
	595	u_p(nzb,-1,:) = u_p(nzb+1,-1,:)
	596	v_p(nzb,0,:) = v_p(nzb+1,0,:)
	597	ENDIF
[1353]	598	w_p(nzb,-1,:) = 0.0_wp
[73]	599
[75]	600	!
[978]	601	!-- Top boundary at the outflow
	602	IF ( ibc_uv_t == 0 ) THEN
	603	u_p(nzt+1,-1,:) = u_init(nzt+1)
	604	v_p(nzt+1,0,:) = v_init(nzt+1)
	605	ELSE
[1742]	606	u_p(nzt+1,-1,:) = u_p(nzt,-1,:)
	607	v_p(nzt+1,0,:) = v_p(nzt,0,:)
[978]	608	ENDIF
[1353]	609	w_p(nzt:nzt+1,-1,:) = 0.0_wp
[978]	610
[75]	611	ENDIF
[73]	612
[75]	613	ENDIF
[73]	614
[3182]	615	IF ( bc_radiation_n ) THEN
[73]	616
[1159]	617	IF ( use_cmax ) THEN
	618	u_p(:,ny+1,:) = u(:,ny,:)
	619	v_p(:,ny+1,:) = v(:,ny,:)
	620	w_p(:,ny+1,:) = w(:,ny,:)
	621	ELSEIF ( .NOT. use_cmax ) THEN
[75]	622
[978]	623	c_max = dy / dt_3d
[75]	624
[1353]	625	c_u_m_l = 0.0_wp
	626	c_v_m_l = 0.0_wp
	627	c_w_m_l = 0.0_wp
[978]	628
[1353]	629	c_u_m = 0.0_wp
	630	c_v_m = 0.0_wp
	631	c_w_m = 0.0_wp
[978]	632
[1]	633	!
[996]	634	!-- Calculate the phase speeds for u, v, and w, first local and then
	635	!-- average along the outflow boundary.
	636	DO k = nzb+1, nzt+1
	637	DO i = nxl, nxr
[73]	638
[106]	639	denom = u_m_n(k,ny,i) - u_m_n(k,ny-1,i)
	640
[1353]	641	IF ( denom /= 0.0_wp ) THEN
[996]	642	c_u(k,i) = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) / ( denom * tsc(2) )
[1353]	643	IF ( c_u(k,i) < 0.0_wp ) THEN
	644	c_u(k,i) = 0.0_wp
[106]	645	ELSEIF ( c_u(k,i) > c_max ) THEN
	646	c_u(k,i) = c_max
	647	ENDIF
	648	ELSE
	649	c_u(k,i) = c_max
[73]	650	ENDIF
	651
[106]	652	denom = v_m_n(k,ny,i) - v_m_n(k,ny-1,i)
[73]	653
[1353]	654	IF ( denom /= 0.0_wp ) THEN
[996]	655	c_v(k,i) = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / ( denom * tsc(2) )
[1353]	656	IF ( c_v(k,i) < 0.0_wp ) THEN
	657	c_v(k,i) = 0.0_wp
[106]	658	ELSEIF ( c_v(k,i) > c_max ) THEN
	659	c_v(k,i) = c_max
	660	ENDIF
	661	ELSE
	662	c_v(k,i) = c_max
[73]	663	ENDIF
	664
[106]	665	denom = w_m_n(k,ny,i) - w_m_n(k,ny-1,i)
[73]	666
[1353]	667	IF ( denom /= 0.0_wp ) THEN
[996]	668	c_w(k,i) = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / ( denom * tsc(2) )
[1353]	669	IF ( c_w(k,i) < 0.0_wp ) THEN
	670	c_w(k,i) = 0.0_wp
[106]	671	ELSEIF ( c_w(k,i) > c_max ) THEN
	672	c_w(k,i) = c_max
	673	ENDIF
	674	ELSE
	675	c_w(k,i) = c_max
[73]	676	ENDIF
[106]	677
[978]	678	c_u_m_l(k) = c_u_m_l(k) + c_u(k,i)
	679	c_v_m_l(k) = c_v_m_l(k) + c_v(k,i)
	680	c_w_m_l(k) = c_w_m_l(k) + c_w(k,i)
[106]	681
[978]	682	ENDDO
	683	ENDDO
[73]	684
[978]	685	#if defined( __parallel )
	686	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
	687	CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
	688	MPI_SUM, comm1dx, ierr )
	689	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
	690	CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
	691	MPI_SUM, comm1dx, ierr )
	692	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
	693	CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
	694	MPI_SUM, comm1dx, ierr )
	695	#else
	696	c_u_m = c_u_m_l
	697	c_v_m = c_v_m_l
	698	c_w_m = c_w_m_l
	699	#endif
	700
	701	c_u_m = c_u_m / (nx+1)
	702	c_v_m = c_v_m / (nx+1)
	703	c_w_m = c_w_m / (nx+1)
	704
[73]	705	!
[978]	706	!-- Save old timelevels for the next timestep
	707	IF ( intermediate_timestep_count == 1 ) THEN
	708	u_m_n(:,:,:) = u(:,ny-1:ny,:)
	709	v_m_n(:,:,:) = v(:,ny-1:ny,:)
	710	w_m_n(:,:,:) = w(:,ny-1:ny,:)
	711	ENDIF
[73]	712
[978]	713	!
	714	!-- Calculate the new velocities
[996]	715	DO k = nzb+1, nzt+1
	716	DO i = nxlg, nxrg
[978]	717	u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * tsc(2) * c_u_m(k) * &
	718	( u(k,ny+1,i) - u(k,ny,i) ) * ddy
[73]	719
[978]	720	v_p(k,ny+1,i) = v(k,ny+1,i) - dt_3d * tsc(2) * c_v_m(k) * &
	721	( v(k,ny+1,i) - v(k,ny,i) ) * ddy
[73]	722
[978]	723	w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * tsc(2) * c_w_m(k) * &
	724	( w(k,ny+1,i) - w(k,ny,i) ) * ddy
	725	ENDDO
[1]	726	ENDDO
	727
	728	!
[978]	729	!-- Bottom boundary at the outflow
	730	IF ( ibc_uv_b == 0 ) THEN
[1353]	731	u_p(nzb,ny+1,:) = 0.0_wp
	732	v_p(nzb,ny+1,:) = 0.0_wp
[978]	733	ELSE
	734	u_p(nzb,ny+1,:) = u_p(nzb+1,ny+1,:)
	735	v_p(nzb,ny+1,:) = v_p(nzb+1,ny+1,:)
	736	ENDIF
[1353]	737	w_p(nzb,ny+1,:) = 0.0_wp
[73]	738
	739	!
[978]	740	!-- Top boundary at the outflow
	741	IF ( ibc_uv_t == 0 ) THEN
	742	u_p(nzt+1,ny+1,:) = u_init(nzt+1)
	743	v_p(nzt+1,ny+1,:) = v_init(nzt+1)
	744	ELSE
	745	u_p(nzt+1,ny+1,:) = u_p(nzt,nyn+1,:)
	746	v_p(nzt+1,ny+1,:) = v_p(nzt,nyn+1,:)
	747	ENDIF
[1353]	748	w_p(nzt:nzt+1,ny+1,:) = 0.0_wp
[978]	749
[1]	750	ENDIF
	751
[75]	752	ENDIF
	753
[3182]	754	IF ( bc_radiation_l ) THEN
[75]	755
[1159]	756	IF ( use_cmax ) THEN
[1717]	757	u_p(:,:,0) = u(:,:,1)
	758	v_p(:,:,-1) = v(:,:,0)
[1159]	759	w_p(:,:,-1) = w(:,:,0)
	760	ELSEIF ( .NOT. use_cmax ) THEN
[75]	761
[978]	762	c_max = dx / dt_3d
[75]	763
[1353]	764	c_u_m_l = 0.0_wp
	765	c_v_m_l = 0.0_wp
	766	c_w_m_l = 0.0_wp
[978]	767
[1353]	768	c_u_m = 0.0_wp
	769	c_v_m = 0.0_wp
	770	c_w_m = 0.0_wp
[978]	771
[1]	772	!
[996]	773	!-- Calculate the phase speeds for u, v, and w, first local and then
	774	!-- average along the outflow boundary.
	775	DO k = nzb+1, nzt+1
	776	DO j = nys, nyn
[75]	777
[106]	778	denom = u_m_l(k,j,1) - u_m_l(k,j,2)
	779
[1353]	780	IF ( denom /= 0.0_wp ) THEN
[996]	781	c_u(k,j) = -c_max * ( u(k,j,1) - u_m_l(k,j,1) ) / ( denom * tsc(2) )
[1353]	782	IF ( c_u(k,j) < 0.0_wp ) THEN
	783	c_u(k,j) = 0.0_wp
[107]	784	ELSEIF ( c_u(k,j) > c_max ) THEN
	785	c_u(k,j) = c_max
[106]	786	ENDIF
	787	ELSE
[107]	788	c_u(k,j) = c_max
[75]	789	ENDIF
	790
[106]	791	denom = v_m_l(k,j,0) - v_m_l(k,j,1)
[75]	792
[1353]	793	IF ( denom /= 0.0_wp ) THEN
[996]	794	c_v(k,j) = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) / ( denom * tsc(2) )
[1353]	795	IF ( c_v(k,j) < 0.0_wp ) THEN
	796	c_v(k,j) = 0.0_wp
[106]	797	ELSEIF ( c_v(k,j) > c_max ) THEN
	798	c_v(k,j) = c_max
	799	ENDIF
	800	ELSE
	801	c_v(k,j) = c_max
[75]	802	ENDIF
	803
[106]	804	denom = w_m_l(k,j,0) - w_m_l(k,j,1)
[75]	805
[1353]	806	IF ( denom /= 0.0_wp ) THEN
[996]	807	c_w(k,j) = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) / ( denom * tsc(2) )
[1353]	808	IF ( c_w(k,j) < 0.0_wp ) THEN
	809	c_w(k,j) = 0.0_wp
[106]	810	ELSEIF ( c_w(k,j) > c_max ) THEN
	811	c_w(k,j) = c_max
	812	ENDIF
	813	ELSE
	814	c_w(k,j) = c_max
[75]	815	ENDIF
[106]	816
[978]	817	c_u_m_l(k) = c_u_m_l(k) + c_u(k,j)
	818	c_v_m_l(k) = c_v_m_l(k) + c_v(k,j)
	819	c_w_m_l(k) = c_w_m_l(k) + c_w(k,j)
[106]	820
[978]	821	ENDDO
	822	ENDDO
[75]	823
[978]	824	#if defined( __parallel )
	825	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
	826	CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
	827	MPI_SUM, comm1dy, ierr )
	828	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
	829	CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
	830	MPI_SUM, comm1dy, ierr )
	831	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
	832	CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
	833	MPI_SUM, comm1dy, ierr )
	834	#else
	835	c_u_m = c_u_m_l
	836	c_v_m = c_v_m_l
	837	c_w_m = c_w_m_l
	838	#endif
	839
	840	c_u_m = c_u_m / (ny+1)
	841	c_v_m = c_v_m / (ny+1)
	842	c_w_m = c_w_m / (ny+1)
	843
[73]	844	!
[978]	845	!-- Save old timelevels for the next timestep
	846	IF ( intermediate_timestep_count == 1 ) THEN
	847	u_m_l(:,:,:) = u(:,:,1:2)
	848	v_m_l(:,:,:) = v(:,:,0:1)
	849	w_m_l(:,:,:) = w(:,:,0:1)
	850	ENDIF
	851
	852	!
	853	!-- Calculate the new velocities
[996]	854	DO k = nzb+1, nzt+1
[1113]	855	DO j = nysg, nyng
[978]	856	u_p(k,j,0) = u(k,j,0) - dt_3d * tsc(2) * c_u_m(k) * &
[106]	857	( u(k,j,0) - u(k,j,1) ) * ddx
[75]	858
[978]	859	v_p(k,j,-1) = v(k,j,-1) - dt_3d * tsc(2) * c_v_m(k) * &
[75]	860	( v(k,j,-1) - v(k,j,0) ) * ddx
	861
[978]	862	w_p(k,j,-1) = w(k,j,-1) - dt_3d * tsc(2) * c_w_m(k) * &
[75]	863	( w(k,j,-1) - w(k,j,0) ) * ddx
[978]	864	ENDDO
[75]	865	ENDDO
	866
	867	!
[978]	868	!-- Bottom boundary at the outflow
	869	IF ( ibc_uv_b == 0 ) THEN
[1353]	870	u_p(nzb,:,0) = 0.0_wp
	871	v_p(nzb,:,-1) = 0.0_wp
[978]	872	ELSE
	873	u_p(nzb,:,0) = u_p(nzb+1,:,0)
	874	v_p(nzb,:,-1) = v_p(nzb+1,:,-1)
	875	ENDIF
[1353]	876	w_p(nzb,:,-1) = 0.0_wp
[1]	877
[75]	878	!
[978]	879	!-- Top boundary at the outflow
	880	IF ( ibc_uv_t == 0 ) THEN
[1764]	881	u_p(nzt+1,:,0) = u_init(nzt+1)
[978]	882	v_p(nzt+1,:,-1) = v_init(nzt+1)
	883	ELSE
[1764]	884	u_p(nzt+1,:,0) = u_p(nzt,:,0)
[978]	885	v_p(nzt+1,:,-1) = v_p(nzt,:,-1)
	886	ENDIF
[1353]	887	w_p(nzt:nzt+1,:,-1) = 0.0_wp
[978]	888
[75]	889	ENDIF
[73]	890
[75]	891	ENDIF
[73]	892
[3182]	893	IF ( bc_radiation_r ) THEN
[73]	894
[1159]	895	IF ( use_cmax ) THEN
	896	u_p(:,:,nx+1) = u(:,:,nx)
	897	v_p(:,:,nx+1) = v(:,:,nx)
	898	w_p(:,:,nx+1) = w(:,:,nx)
	899	ELSEIF ( .NOT. use_cmax ) THEN
[75]	900
[978]	901	c_max = dx / dt_3d
[75]	902
[1353]	903	c_u_m_l = 0.0_wp
	904	c_v_m_l = 0.0_wp
	905	c_w_m_l = 0.0_wp
[978]	906
[1353]	907	c_u_m = 0.0_wp
	908	c_v_m = 0.0_wp
	909	c_w_m = 0.0_wp
[978]	910
[1]	911	!
[996]	912	!-- Calculate the phase speeds for u, v, and w, first local and then
	913	!-- average along the outflow boundary.
	914	DO k = nzb+1, nzt+1
	915	DO j = nys, nyn
[73]	916
[106]	917	denom = u_m_r(k,j,nx) - u_m_r(k,j,nx-1)
	918
[1353]	919	IF ( denom /= 0.0_wp ) THEN
[996]	920	c_u(k,j) = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) ) / ( denom * tsc(2) )
[1353]	921	IF ( c_u(k,j) < 0.0_wp ) THEN
	922	c_u(k,j) = 0.0_wp
[106]	923	ELSEIF ( c_u(k,j) > c_max ) THEN
	924	c_u(k,j) = c_max
	925	ENDIF
	926	ELSE
	927	c_u(k,j) = c_max
[73]	928	ENDIF
	929
[106]	930	denom = v_m_r(k,j,nx) - v_m_r(k,j,nx-1)
[73]	931
[1353]	932	IF ( denom /= 0.0_wp ) THEN
[996]	933	c_v(k,j) = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) / ( denom * tsc(2) )
[1353]	934	IF ( c_v(k,j) < 0.0_wp ) THEN
	935	c_v(k,j) = 0.0_wp
[106]	936	ELSEIF ( c_v(k,j) > c_max ) THEN
	937	c_v(k,j) = c_max
	938	ENDIF
	939	ELSE
	940	c_v(k,j) = c_max
[73]	941	ENDIF
	942
[106]	943	denom = w_m_r(k,j,nx) - w_m_r(k,j,nx-1)
[73]	944
[1353]	945	IF ( denom /= 0.0_wp ) THEN
[996]	946	c_w(k,j) = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) / ( denom * tsc(2) )
[1353]	947	IF ( c_w(k,j) < 0.0_wp ) THEN
	948	c_w(k,j) = 0.0_wp
[106]	949	ELSEIF ( c_w(k,j) > c_max ) THEN
	950	c_w(k,j) = c_max
	951	ENDIF
	952	ELSE
	953	c_w(k,j) = c_max
[73]	954	ENDIF
[106]	955
[978]	956	c_u_m_l(k) = c_u_m_l(k) + c_u(k,j)
	957	c_v_m_l(k) = c_v_m_l(k) + c_v(k,j)
	958	c_w_m_l(k) = c_w_m_l(k) + c_w(k,j)
[106]	959
[978]	960	ENDDO
	961	ENDDO
[73]	962
[978]	963	#if defined( __parallel )
	964	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
	965	CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
	966	MPI_SUM, comm1dy, ierr )
	967	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
	968	CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
	969	MPI_SUM, comm1dy, ierr )
	970	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
	971	CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
	972	MPI_SUM, comm1dy, ierr )
	973	#else
	974	c_u_m = c_u_m_l
	975	c_v_m = c_v_m_l
	976	c_w_m = c_w_m_l
	977	#endif
	978
	979	c_u_m = c_u_m / (ny+1)
	980	c_v_m = c_v_m / (ny+1)
	981	c_w_m = c_w_m / (ny+1)
	982
[73]	983	!
[978]	984	!-- Save old timelevels for the next timestep
	985	IF ( intermediate_timestep_count == 1 ) THEN
	986	u_m_r(:,:,:) = u(:,:,nx-1:nx)
	987	v_m_r(:,:,:) = v(:,:,nx-1:nx)
	988	w_m_r(:,:,:) = w(:,:,nx-1:nx)
	989	ENDIF
[73]	990
[978]	991	!
	992	!-- Calculate the new velocities
[996]	993	DO k = nzb+1, nzt+1
[1113]	994	DO j = nysg, nyng
[978]	995	u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * tsc(2) * c_u_m(k) * &
	996	( u(k,j,nx+1) - u(k,j,nx) ) * ddx
[73]	997
[978]	998	v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * tsc(2) * c_v_m(k) * &
	999	( v(k,j,nx+1) - v(k,j,nx) ) * ddx
[73]	1000
[978]	1001	w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * tsc(2) * c_w_m(k) * &
	1002	( w(k,j,nx+1) - w(k,j,nx) ) * ddx
	1003	ENDDO
[73]	1004	ENDDO
	1005
	1006	!
[978]	1007	!-- Bottom boundary at the outflow
	1008	IF ( ibc_uv_b == 0 ) THEN
[1353]	1009	u_p(nzb,:,nx+1) = 0.0_wp
	1010	v_p(nzb,:,nx+1) = 0.0_wp
[978]	1011	ELSE
	1012	u_p(nzb,:,nx+1) = u_p(nzb+1,:,nx+1)
	1013	v_p(nzb,:,nx+1) = v_p(nzb+1,:,nx+1)
	1014	ENDIF
[1353]	1015	w_p(nzb,:,nx+1) = 0.0_wp
[73]	1016
	1017	!
[978]	1018	!-- Top boundary at the outflow
	1019	IF ( ibc_uv_t == 0 ) THEN
	1020	u_p(nzt+1,:,nx+1) = u_init(nzt+1)
	1021	v_p(nzt+1,:,nx+1) = v_init(nzt+1)
	1022	ELSE
	1023	u_p(nzt+1,:,nx+1) = u_p(nzt,:,nx+1)
	1024	v_p(nzt+1,:,nx+1) = v_p(nzt,:,nx+1)
	1025	ENDIF
[1742]	1026	w_p(nzt:nzt+1,:,nx+1) = 0.0_wp
[978]	1027
[1]	1028	ENDIF
	1029
	1030	ENDIF
[3864]	1031
[3467]	1032	IF ( salsa ) THEN
	1033	CALL salsa_boundary_conds
[3864]	1034	ENDIF
[1]	1035
	1036	END SUBROUTINE boundary_conds

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |