Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |