Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |