Home

Context Navigation

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

Last change on this file since 107 was 107, checked in by raasch, 17 years ago
further bugfix for non-cyclic BCs
Property svn:keywords set to `Id`
File size: 19.5 KB

Rev	Line
[1]	1	SUBROUTINE boundary_conds( range )
	2
	3	!------------------------------------------------------------------------------!
	4	! Actual revisions:
	5	! -----------------
[106]	6	! Boundary conditions for temperature adjusted for coupled runs,
	7	! bugfixes for the radiation boundary conditions at the outflow: radiation
	8	! conditions are used for every substep, phase speeds are calculated for the
	9	! first Runge-Kutta substep only and then reused, several index values changed
	10	!
[1]	11	!
	12	! Former revisions:
	13	! -----------------
[3]	14	! $Id: boundary_conds.f90 107 2007-08-17 13:54:45Z raasch $
[39]	15	!
[98]	16	! 95 2007-06-02 16:48:38Z raasch
	17	! Boundary conditions for salinity added
	18	!
[77]	19	! 75 2007-03-22 09:54:05Z raasch
	20	! The "main" part sets conditions for time level t+dt instead of level t,
	21	! outflow boundary conditions changed from Neumann to radiation condition,
	22	! uxrp, vynp eliminated, moisture renamed humidity
	23	!
[39]	24	! 19 2007-02-23 04:53:48Z raasch
	25	! Boundary conditions for e(nzt), pt(nzt), and q(nzt) removed because these
	26	! gridpoints are now calculated by the prognostic equation,
	27	! Dirichlet and zero gradient condition for pt established at top boundary
	28	!
[3]	29	! RCS Log replace by Id keyword, revision history cleaned up
	30	!
[1]	31	! Revision 1.15 2006/02/23 09:54:55 raasch
	32	! Surface boundary conditions in case of topography: nzb replaced by
	33	! 2d-k-index-arrays (nzb_w_inner, etc.). Conditions for u and v remain
	34	! unchanged (still using nzb) because a non-flat topography must use a
	35	! Prandtl-layer, which don't requires explicit setting of the surface values.
	36	!
	37	! Revision 1.1 1997/09/12 06:21:34 raasch
	38	! Initial revision
	39	!
	40	!
	41	! Description:
	42	! ------------
	43	! Boundary conditions for the prognostic quantities (range='main').
	44	! In case of non-cyclic lateral boundaries the conditions for velocities at
	45	! the outflow are set after the pressure solver has been called (range=
	46	! 'outflow_uvw').
	47	! One additional bottom boundary condition is applied for the TKE (=(u)*2)
	48	! in prandtl_fluxes. The cyclic lateral boundary conditions are implicitly
	49	! handled in routine exchange_horiz. Pressure boundary conditions are
	50	! explicitly set in routines pres, poisfft, poismg and sor.
	51	!------------------------------------------------------------------------------!
	52
	53	USE arrays_3d
	54	USE control_parameters
	55	USE grid_variables
	56	USE indices
	57	USE pegrid
	58
	59	IMPLICIT NONE
	60
	61	CHARACTER (LEN=*) :: range
	62
	63	INTEGER :: i, j, k
	64
[106]	65	REAL :: c_max, denom
[1]	66
[73]	67
[1]	68	IF ( range == 'main') THEN
	69	!
	70	!-- Bottom boundary
	71	IF ( ibc_uv_b == 0 ) THEN
[73]	72	!
	73	!-- Satisfying the Dirichlet condition with an extra layer below the
	74	!-- surface where the u and v component change their sign
	75	u_p(nzb,:,:) = -u_p(nzb+1,:,:)
	76	v_p(nzb,:,:) = -v_p(nzb+1,:,:)
	77	ELSE
	78	u_p(nzb,:,:) = u_p(nzb+1,:,:)
	79	v_p(nzb,:,:) = v_p(nzb+1,:,:)
[1]	80	ENDIF
	81	DO i = nxl-1, nxr+1
	82	DO j = nys-1, nyn+1
[73]	83	w_p(nzb_w_inner(j,i),j,i) = 0.0
[1]	84	ENDDO
	85	ENDDO
	86
	87	!
	88	!-- Top boundary
	89	IF ( ibc_uv_t == 0 ) THEN
[73]	90	u_p(nzt+1,:,:) = ug(nzt+1)
	91	v_p(nzt+1,:,:) = vg(nzt+1)
[1]	92	ELSE
[73]	93	u_p(nzt+1,:,:) = u_p(nzt,:,:)
	94	v_p(nzt+1,:,:) = v_p(nzt,:,:)
[1]	95	ENDIF
[73]	96	w_p(nzt:nzt+1,:,:) = 0.0 ! nzt is not a prognostic level (but cf. pres)
[1]	97
	98	!
[102]	99	!-- Temperature at bottom boundary.
	100	!-- In case of coupled runs (ibc_pt_b = 2) the temperature is given by
	101	!-- the sea surface temperature of the coupled ocean model.
[1]	102	IF ( ibc_pt_b == 0 ) THEN
[73]	103	DO i = nxl-1, nxr+1
	104	DO j = nys-1, nyn+1
	105	pt_p(nzb_s_inner(j,i),j,i) = pt(nzb_s_inner(j,i),j,i)
[1]	106	ENDDO
[73]	107	ENDDO
[102]	108	ELSEIF ( ibc_pt_b == 1 ) THEN
[1]	109	DO i = nxl-1, nxr+1
	110	DO j = nys-1, nyn+1
[73]	111	pt_p(nzb_s_inner(j,i),j,i) = pt_p(nzb_s_inner(j,i)+1,j,i)
[1]	112	ENDDO
	113	ENDDO
	114	ENDIF
	115
	116	!
	117	!-- Temperature at top boundary
[19]	118	IF ( ibc_pt_t == 0 ) THEN
[73]	119	pt_p(nzt+1,:,:) = pt(nzt+1,:,:)
[19]	120	ELSEIF ( ibc_pt_t == 1 ) THEN
[73]	121	pt_p(nzt+1,:,:) = pt_p(nzt,:,:)
[19]	122	ELSEIF ( ibc_pt_t == 2 ) THEN
[73]	123	pt_p(nzt+1,:,:) = pt_p(nzt,:,:) + bc_pt_t_val * dzu(nzt+1)
[1]	124	ENDIF
	125
	126	!
	127	!-- Boundary conditions for TKE
	128	!-- Generally Neumann conditions with de/dz=0 are assumed
	129	IF ( .NOT. constant_diffusion ) THEN
	130	DO i = nxl-1, nxr+1
	131	DO j = nys-1, nyn+1
[73]	132	e_p(nzb_s_inner(j,i),j,i) = e_p(nzb_s_inner(j,i)+1,j,i)
[1]	133	ENDDO
	134	ENDDO
[73]	135	e_p(nzt+1,:,:) = e_p(nzt,:,:)
[1]	136	ENDIF
	137
	138	!
[95]	139	!-- Boundary conditions for salinity
	140	IF ( ocean ) THEN
	141	!
	142	!-- Bottom boundary: Neumann condition because salinity flux is always
	143	!-- given
	144	DO i = nxl-1, nxr+1
	145	DO j = nys-1, nyn+1
	146	sa_p(nzb_s_inner(j,i),j,i) = sa_p(nzb_s_inner(j,i)+1,j,i)
	147	ENDDO
	148	ENDDO
	149
	150	!
	151	!-- Top boundary: Dirichlet or Neumann
	152	IF ( ibc_sa_t == 0 ) THEN
	153	sa_p(nzt+1,:,:) = sa(nzt+1,:,:)
	154	ELSEIF ( ibc_sa_t == 1 ) THEN
	155	sa_p(nzt+1,:,:) = sa_p(nzt,:,:)
	156	ENDIF
	157
	158	ENDIF
	159
	160	!
[1]	161	!-- Boundary conditions for total water content or scalar,
[95]	162	!-- bottom and top boundary (see also temperature)
[75]	163	IF ( humidity .OR. passive_scalar ) THEN
[1]	164	!
[75]	165	!-- Surface conditions for constant_humidity_flux
[1]	166	IF ( ibc_q_b == 0 ) THEN
[73]	167	DO i = nxl-1, nxr+1
	168	DO j = nys-1, nyn+1
	169	q_p(nzb_s_inner(j,i),j,i) = q(nzb_s_inner(j,i),j,i)
[1]	170	ENDDO
[73]	171	ENDDO
[1]	172	ELSE
	173	DO i = nxl-1, nxr+1
	174	DO j = nys-1, nyn+1
[73]	175	q_p(nzb_s_inner(j,i),j,i) = q_p(nzb_s_inner(j,i)+1,j,i)
[1]	176	ENDDO
	177	ENDDO
	178	ENDIF
	179	!
	180	!-- Top boundary
[73]	181	q_p(nzt+1,:,:) = q_p(nzt,:,:) + bc_q_t_val * dzu(nzt+1)
[1]	182	ENDIF
	183
	184	!
	185	!-- Lateral boundary conditions at the inflow. Quasi Neumann conditions
	186	!-- are needed for the wall normal velocity in order to ensure zero
	187	!-- divergence. Dirichlet conditions are used for all other quantities.
	188	IF ( inflow_s ) THEN
[73]	189	v_p(:,nys,:) = v_p(:,nys-1,:)
[1]	190	ELSEIF ( inflow_n ) THEN
[75]	191	v_p(:,nyn,:) = v_p(:,nyn+1,:)
[1]	192	ELSEIF ( inflow_l ) THEN
[73]	193	u_p(:,:,nxl) = u_p(:,:,nxl-1)
[1]	194	ELSEIF ( inflow_r ) THEN
[75]	195	u_p(:,:,nxr) = u_p(:,:,nxr+1)
[1]	196	ENDIF
	197
	198	!
	199	!-- Lateral boundary conditions for scalar quantities at the outflow
	200	IF ( outflow_s ) THEN
[73]	201	pt_p(:,nys-1,:) = pt_p(:,nys,:)
	202	IF ( .NOT. constant_diffusion ) e_p(:,nys-1,:) = e_p(:,nys,:)
[75]	203	IF ( humidity .OR. passive_scalar ) q_p(:,nys-1,:) = q_p(:,nys,:)
[1]	204	ELSEIF ( outflow_n ) THEN
[73]	205	pt_p(:,nyn+1,:) = pt_p(:,nyn,:)
	206	IF ( .NOT. constant_diffusion ) e_p(:,nyn+1,:) = e_p(:,nyn,:)
[75]	207	IF ( humidity .OR. passive_scalar ) q_p(:,nyn+1,:) = q_p(:,nyn,:)
[1]	208	ELSEIF ( outflow_l ) THEN
[73]	209	pt_p(:,:,nxl-1) = pt_p(:,:,nxl)
	210	IF ( .NOT. constant_diffusion ) e_p(:,:,nxl-1) = e_p(:,:,nxl)
[75]	211	IF ( humidity .OR. passive_scalar ) q_p(:,:,nxl-1) = q_p(:,:,nxl)
[1]	212	ELSEIF ( outflow_r ) THEN
[73]	213	pt_p(:,:,nxr+1) = pt_p(:,:,nxr)
	214	IF ( .NOT. constant_diffusion ) e_p(:,:,nxr+1) = e_p(:,:,nxr)
[75]	215	IF ( humidity .OR. passive_scalar ) q_p(:,:,nxr+1) = q_p(:,:,nxr)
[1]	216	ENDIF
	217
	218	ENDIF
	219
	220	!
[75]	221	!-- Radiation boundary condition for the velocities at the respective outflow
[106]	222	IF ( outflow_s ) THEN
[75]	223
	224	c_max = dy / dt_3d
	225
	226	DO i = nxl-1, nxr+1
	227	DO k = nzb+1, nzt+1
	228
	229	!
[106]	230	!-- Calculate the phase speeds for u,v, and w. In case of using a
	231	!-- Runge-Kutta scheme, do this for the first substep only and then
	232	!-- reuse this values for the further substeps.
	233	IF ( intermediate_timestep_count == 1 ) THEN
[75]	234
[106]	235	denom = u_m_s(k,0,i) - u_m_s(k,1,i)
	236
	237	IF ( denom /= 0.0 ) THEN
	238	c_u(k,i) = -c_max * ( u(k,0,i) - u_m_s(k,0,i) ) / denom
	239	IF ( c_u(k,i) < 0.0 ) THEN
	240	c_u(k,i) = 0.0
	241	ELSEIF ( c_u(k,i) > c_max ) THEN
	242	c_u(k,i) = c_max
	243	ENDIF
	244	ELSE
	245	c_u(k,i) = c_max
[75]	246	ENDIF
	247
[106]	248	denom = v_m_s(k,1,i) - v_m_s(k,2,i)
	249
	250	IF ( denom /= 0.0 ) THEN
	251	c_v(k,i) = -c_max * ( v(k,1,i) - v_m_s(k,1,i) ) / denom
	252	IF ( c_v(k,i) < 0.0 ) THEN
	253	c_v(k,i) = 0.0
	254	ELSEIF ( c_v(k,i) > c_max ) THEN
	255	c_v(k,i) = c_max
	256	ENDIF
	257	ELSE
	258	c_v(k,i) = c_max
[75]	259	ENDIF
	260
[106]	261	denom = w_m_s(k,0,i) - w_m_s(k,1,i)
[75]	262
[106]	263	IF ( denom /= 0.0 ) THEN
	264	c_w(k,i) = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) / denom
	265	IF ( c_w(k,i) < 0.0 ) THEN
	266	c_w(k,i) = 0.0
	267	ELSEIF ( c_w(k,i) > c_max ) THEN
	268	c_w(k,i) = c_max
	269	ENDIF
	270	ELSE
	271	c_w(k,i) = c_max
[75]	272	ENDIF
[106]	273
	274	!
	275	!-- Save old timelevels for the next timestep
	276	u_m_s(k,:,i) = u(k,0:1,i)
	277	v_m_s(k,:,i) = v(k,1:2,i)
	278	w_m_s(k,:,i) = w(k,0:1,i)
	279
[75]	280	ENDIF
	281
	282	!
	283	!-- Calculate the new velocities
[106]	284	u_p(k,-1,i) = u(k,-1,i) - dt_3d * tsc(2) * c_u(k,i) * &
[75]	285	( u(k,-1,i) - u(k,0,i) ) * ddy
	286
[107]	287	v_p(k,0,i) = v(k,0,i) - dt_3d * tsc(2) * c_v(k,i) * &
[106]	288	( v(k,0,i) - v(k,1,i) ) * ddy
[75]	289
[106]	290	w_p(k,-1,i) = w(k,-1,i) - dt_3d * tsc(2) * c_w(k,i) * &
[75]	291	( w(k,-1,i) - w(k,0,i) ) * ddy
	292
	293	ENDDO
	294	ENDDO
	295
	296	!
	297	!-- Bottom boundary at the outflow
	298	IF ( ibc_uv_b == 0 ) THEN
	299	u_p(nzb,-1,:) = -u_p(nzb+1,-1,:)
[106]	300	v_p(nzb,0,:) = -v_p(nzb+1,0,:)
[75]	301	ELSE
	302	u_p(nzb,-1,:) = u_p(nzb+1,-1,:)
[106]	303	v_p(nzb,0,:) = v_p(nzb+1,0,:)
[73]	304	ENDIF
[106]	305	w_p(nzb,-1,:) = 0.0
[73]	306
[75]	307	!
	308	!-- Top boundary at the outflow
	309	IF ( ibc_uv_t == 0 ) THEN
	310	u_p(nzt+1,-1,:) = ug(nzt+1)
[106]	311	v_p(nzt+1,0,:) = vg(nzt+1)
[75]	312	ELSE
	313	u_p(nzt+1,-1,:) = u(nzt,-1,:)
[106]	314	v_p(nzt+1,0,:) = v(nzt,0,:)
[75]	315	ENDIF
	316	w_p(nzt:nzt+1,-1,:) = 0.0
[73]	317
[75]	318	ENDIF
[73]	319
[106]	320	IF ( outflow_n ) THEN
[73]	321
[75]	322	c_max = dy / dt_3d
	323
	324	DO i = nxl-1, nxr+1
	325	DO k = nzb+1, nzt+1
	326
[1]	327	!
[106]	328	!-- Calculate the phase speeds for u,v, and w. In case of using a
	329	!-- Runge-Kutta scheme, do this for the first substep only and then
	330	!-- reuse this values for the further substeps.
	331	IF ( intermediate_timestep_count == 1 ) THEN
[73]	332
[106]	333	denom = u_m_n(k,ny,i) - u_m_n(k,ny-1,i)
	334
	335	IF ( denom /= 0.0 ) THEN
	336	c_u(k,i) = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) / denom
	337	IF ( c_u(k,i) < 0.0 ) THEN
	338	c_u(k,i) = 0.0
	339	ELSEIF ( c_u(k,i) > c_max ) THEN
	340	c_u(k,i) = c_max
	341	ENDIF
	342	ELSE
	343	c_u(k,i) = c_max
[73]	344	ENDIF
	345
[106]	346	denom = v_m_n(k,ny,i) - v_m_n(k,ny-1,i)
[73]	347
[106]	348	IF ( denom /= 0.0 ) THEN
	349	c_v(k,i) = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / denom
	350	IF ( c_v(k,i) < 0.0 ) THEN
	351	c_v(k,i) = 0.0
	352	ELSEIF ( c_v(k,i) > c_max ) THEN
	353	c_v(k,i) = c_max
	354	ENDIF
	355	ELSE
	356	c_v(k,i) = c_max
[73]	357	ENDIF
	358
[106]	359	denom = w_m_n(k,ny,i) - w_m_n(k,ny-1,i)
[73]	360
[106]	361	IF ( denom /= 0.0 ) THEN
	362	c_w(k,i) = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / denom
	363	IF ( c_w(k,i) < 0.0 ) THEN
	364	c_w(k,i) = 0.0
	365	ELSEIF ( c_w(k,i) > c_max ) THEN
	366	c_w(k,i) = c_max
	367	ENDIF
	368	ELSE
	369	c_w(k,i) = c_max
[73]	370	ENDIF
[106]	371
	372	!
	373	!-- Swap timelevels for the next timestep
	374	u_m_n(k,:,i) = u(k,ny-1:ny,i)
	375	v_m_n(k,:,i) = v(k,ny-1:ny,i)
	376	w_m_n(k,:,i) = w(k,ny-1:ny,i)
	377
[75]	378	ENDIF
[73]	379
	380	!
[75]	381	!-- Calculate the new velocities
[106]	382	u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * tsc(2) * c_u(k,i) * &
[75]	383	( u(k,ny+1,i) - u(k,ny,i) ) * ddy
[73]	384
[106]	385	v_p(k,ny+1,i) = v(k,ny+1,i) - dt_3d * tsc(2) * c_v(k,i) * &
[75]	386	( v(k,ny+1,i) - v(k,ny,i) ) * ddy
[73]	387
[106]	388	w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * tsc(2) * c_w(k,i) * &
[75]	389	( w(k,ny+1,i) - w(k,ny,i) ) * ddy
[73]	390
[1]	391	ENDDO
[75]	392	ENDDO
[1]	393
	394	!
[75]	395	!-- Bottom boundary at the outflow
	396	IF ( ibc_uv_b == 0 ) THEN
	397	u_p(nzb,ny+1,:) = -u_p(nzb+1,ny+1,:)
	398	v_p(nzb,ny+1,:) = -v_p(nzb+1,ny+1,:)
	399	ELSE
	400	u_p(nzb,ny+1,:) = u_p(nzb+1,ny+1,:)
	401	v_p(nzb,ny+1,:) = v_p(nzb+1,ny+1,:)
	402	ENDIF
	403	w_p(nzb,ny+1,:) = 0.0
[73]	404
	405	!
[75]	406	!-- Top boundary at the outflow
	407	IF ( ibc_uv_t == 0 ) THEN
	408	u_p(nzt+1,ny+1,:) = ug(nzt+1)
	409	v_p(nzt+1,ny+1,:) = vg(nzt+1)
	410	ELSE
	411	u_p(nzt+1,ny+1,:) = u_p(nzt,nyn+1,:)
	412	v_p(nzt+1,ny+1,:) = v_p(nzt,nyn+1,:)
[1]	413	ENDIF
[75]	414	w_p(nzt:nzt+1,ny+1,:) = 0.0
[1]	415
[75]	416	ENDIF
	417
[106]	418	IF ( outflow_l ) THEN
[75]	419
	420	c_max = dx / dt_3d
	421
	422	DO j = nys-1, nyn+1
	423	DO k = nzb+1, nzt+1
	424
[1]	425	!
[106]	426	!-- Calculate the phase speeds for u,v, and w. In case of using a
	427	!-- Runge-Kutta scheme, do this for the first substep only and then
	428	!-- reuse this values for the further substeps.
	429	IF ( intermediate_timestep_count == 1 ) THEN
[75]	430
[106]	431	denom = u_m_l(k,j,1) - u_m_l(k,j,2)
	432
	433	IF ( denom /= 0.0 ) THEN
	434	c_u(k,j) = -c_max * ( u(k,j,1) - u_m_l(k,j,1) ) / denom
[107]	435	IF ( c_u(k,j) < 0.0 ) THEN
[106]	436	c_u(k,j) = 0.0
[107]	437	ELSEIF ( c_u(k,j) > c_max ) THEN
	438	c_u(k,j) = c_max
[106]	439	ENDIF
	440	ELSE
[107]	441	c_u(k,j) = c_max
[75]	442	ENDIF
	443
[106]	444	denom = v_m_l(k,j,0) - v_m_l(k,j,1)
[75]	445
[106]	446	IF ( denom /= 0.0 ) THEN
	447	c_v(k,j) = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) / denom
	448	IF ( c_v(k,j) < 0.0 ) THEN
	449	c_v(k,j) = 0.0
	450	ELSEIF ( c_v(k,j) > c_max ) THEN
	451	c_v(k,j) = c_max
	452	ENDIF
	453	ELSE
	454	c_v(k,j) = c_max
[75]	455	ENDIF
	456
[106]	457	denom = w_m_l(k,j,0) - w_m_l(k,j,1)
[75]	458
[106]	459	IF ( denom /= 0.0 ) THEN
	460	c_w(k,j) = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) / denom
	461	IF ( c_w(k,j) < 0.0 ) THEN
	462	c_w(k,j) = 0.0
	463	ELSEIF ( c_w(k,j) > c_max ) THEN
	464	c_w(k,j) = c_max
	465	ENDIF
	466	ELSE
	467	c_w(k,j) = c_max
[75]	468	ENDIF
[106]	469
	470	!
	471	!-- Swap timelevels for the next timestep
	472	u_m_l(k,j,:) = u(k,j,1:2)
	473	v_m_l(k,j,:) = v(k,j,0:1)
	474	w_m_l(k,j,:) = w(k,j,0:1)
	475
[75]	476	ENDIF
	477
[73]	478	!
[75]	479	!-- Calculate the new velocities
[106]	480	u_p(k,j,0) = u(k,j,0) - dt_3d * tsc(2) * c_u(k,j) * &
	481	( u(k,j,0) - u(k,j,1) ) * ddx
[75]	482
[106]	483	v_p(k,j,-1) = v(k,j,-1) - dt_3d * tsc(2) * c_v(k,j) * &
[75]	484	( v(k,j,-1) - v(k,j,0) ) * ddx
	485
[106]	486	w_p(k,j,-1) = w(k,j,-1) - dt_3d * tsc(2) * c_w(k,j) * &
[75]	487	( w(k,j,-1) - w(k,j,0) ) * ddx
	488
	489	ENDDO
	490	ENDDO
	491
	492	!
	493	!-- Bottom boundary at the outflow
	494	IF ( ibc_uv_b == 0 ) THEN
	495	u_p(nzb,:,-1) = -u_p(nzb+1,:,-1)
	496	v_p(nzb,:,-1) = -v_p(nzb+1,:,-1)
	497	ELSE
	498	u_p(nzb,:,-1) = u_p(nzb+1,:,-1)
	499	v_p(nzb,:,-1) = v_p(nzb+1,:,-1)
[1]	500	ENDIF
[75]	501	w_p(nzb,:,-1) = 0.0
[1]	502
[75]	503	!
	504	!-- Top boundary at the outflow
	505	IF ( ibc_uv_t == 0 ) THEN
	506	u_p(nzt+1,:,-1) = ug(nzt+1)
	507	v_p(nzt+1,:,-1) = vg(nzt+1)
	508	ELSE
	509	u_p(nzt+1,:,-1) = u_p(nzt,:,-1)
	510	v_p(nzt+1,:,-1) = v_p(nzt,:,-1)
	511	ENDIF
	512	w_p(nzt:nzt+1,:,-1) = 0.0
[73]	513
[75]	514	ENDIF
[73]	515
[106]	516	IF ( outflow_r ) THEN
[73]	517
[75]	518	c_max = dx / dt_3d
	519
	520	DO j = nys-1, nyn+1
	521	DO k = nzb+1, nzt+1
	522
[1]	523	!
[106]	524	!-- Calculate the phase speeds for u,v, and w. In case of using a
	525	!-- Runge-Kutta scheme, do this for the first substep only and then
	526	!-- reuse this values for the further substeps.
	527	IF ( intermediate_timestep_count == 1 ) THEN
[73]	528
[106]	529	denom = u_m_r(k,j,nx) - u_m_r(k,j,nx-1)
	530
	531	IF ( denom /= 0.0 ) THEN
	532	c_u(k,j) = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) ) / denom
	533	IF ( c_u(k,j) < 0.0 ) THEN
	534	c_u(k,j) = 0.0
	535	ELSEIF ( c_u(k,j) > c_max ) THEN
	536	c_u(k,j) = c_max
	537	ENDIF
	538	ELSE
	539	c_u(k,j) = c_max
[73]	540	ENDIF
	541
[106]	542	denom = v_m_r(k,j,nx) - v_m_r(k,j,nx-1)
[73]	543
[106]	544	IF ( denom /= 0.0 ) THEN
	545	c_v(k,j) = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) / denom
	546	IF ( c_v(k,j) < 0.0 ) THEN
	547	c_v(k,j) = 0.0
	548	ELSEIF ( c_v(k,j) > c_max ) THEN
	549	c_v(k,j) = c_max
	550	ENDIF
	551	ELSE
	552	c_v(k,j) = c_max
[73]	553	ENDIF
	554
[106]	555	denom = w_m_r(k,j,nx) - w_m_r(k,j,nx-1)
[73]	556
[106]	557	IF ( denom /= 0.0 ) THEN
	558	c_w(k,j) = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) / denom
	559	IF ( c_w(k,j) < 0.0 ) THEN
	560	c_w(k,j) = 0.0
	561	ELSEIF ( c_w(k,j) > c_max ) THEN
	562	c_w(k,j) = c_max
	563	ENDIF
	564	ELSE
	565	c_w(k,j) = c_max
[73]	566	ENDIF
[106]	567
	568	!
	569	!-- Swap timelevels for the next timestep
	570	u_m_r(k,j,:) = u(k,j,nx-1:nx)
	571	v_m_r(k,j,:) = v(k,j,nx-1:nx)
	572	w_m_r(k,j,:) = w(k,j,nx-1:nx)
	573
[75]	574	ENDIF
[73]	575
	576	!
[75]	577	!-- Calculate the new velocities
[106]	578	u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * tsc(2) * c_u(k,j) * &
[75]	579	( u(k,j,nx+1) - u(k,j,nx) ) * ddx
[73]	580
[106]	581	v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * tsc(2) * c_v(k,j) * &
[75]	582	( v(k,j,nx+1) - v(k,j,nx) ) * ddx
[73]	583
[106]	584	w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * tsc(2) * c_w(k,j) * &
[75]	585	( w(k,j,nx+1) - w(k,j,nx) ) * ddx
[73]	586
	587	ENDDO
[75]	588	ENDDO
[73]	589
	590	!
[75]	591	!-- Bottom boundary at the outflow
	592	IF ( ibc_uv_b == 0 ) THEN
	593	u_p(nzb,:,nx+1) = -u_p(nzb+1,:,nx+1)
	594	v_p(nzb,:,nx+1) = -v_p(nzb+1,:,nx+1)
	595	ELSE
	596	u_p(nzb,:,nx+1) = u_p(nzb+1,:,nx+1)
	597	v_p(nzb,:,nx+1) = v_p(nzb+1,:,nx+1)
	598	ENDIF
	599	w_p(nzb,:,nx+1) = 0.0
[73]	600
	601	!
[75]	602	!-- Top boundary at the outflow
	603	IF ( ibc_uv_t == 0 ) THEN
	604	u_p(nzt+1,:,nx+1) = ug(nzt+1)
	605	v_p(nzt+1,:,nx+1) = vg(nzt+1)
	606	ELSE
	607	u_p(nzt+1,:,nx+1) = u_p(nzt,:,nx+1)
	608	v_p(nzt+1,:,nx+1) = v_p(nzt,:,nx+1)
[1]	609	ENDIF
[75]	610	w(nzt:nzt+1,:,nx+1) = 0.0
[1]	611
	612	ENDIF
	613
	614
	615	END SUBROUTINE boundary_conds

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |