Home

Context Navigation

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

Last change on this file since 95 was 95, checked in by raasch, 17 years ago
further preliminary uncomplete changes for ocean version
Property svn:keywords set to `Id`
File size: 17.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |