Home

Context Navigation

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

Last change on this file since 1025 was 997, checked in by raasch, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 26.1 KB

Rev	Line
[1]	1	SUBROUTINE boundary_conds( range )
	2
	3	!------------------------------------------------------------------------------!
[484]	4	! Current revisions:
[1]	5	! -----------------
[997]	6	!
[667]	7	!
[1]	8	! Former revisions:
	9	! -----------------
[3]	10	! $Id: boundary_conds.f90 997 2012-09-07 10:53:03Z letzel $
[39]	11	!
[997]	12	! 996 2012-09-07 10:41:47Z raasch
	13	! little reformatting
	14	!
[979]	15	! 978 2012-08-09 08:28:32Z fricke
	16	! Neumann boudnary conditions are added at the inflow boundary for the SGS-TKE.
	17	! Outflow boundary conditions for the velocity components can be set to Neumann
	18	! conditions or to radiation conditions with a horizontal averaged phase
	19	! velocity.
	20	!
[876]	21	! 875 2012-04-02 15:35:15Z gryschka
	22	! Bugfix in case of dirichlet inflow bc at the right or north boundary
	23	!
[768]	24	! 767 2011-10-14 06:39:12Z raasch
	25	! ug,vg replaced by u_init,v_init as the Dirichlet top boundary condition
	26	!
[668]	27	! 667 2010-12-23 12:06:00Z suehring/gryschka
	28	! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
	29	! Removed mirror boundary conditions for u and v at the bottom in case of
	30	! ibc_uv_b == 0. Instead, dirichelt boundary conditions (u=v=0) are set
	31	! in init_3d_model
	32	!
[110]	33	! 107 2007-08-17 13:54:45Z raasch
	34	! Boundary conditions for temperature adjusted for coupled runs,
	35	! bugfixes for the radiation boundary conditions at the outflow: radiation
	36	! conditions are used for every substep, phase speeds are calculated for the
	37	! first Runge-Kutta substep only and then reused, several index values changed
	38	!
[98]	39	! 95 2007-06-02 16:48:38Z raasch
	40	! Boundary conditions for salinity added
	41	!
[77]	42	! 75 2007-03-22 09:54:05Z raasch
	43	! The "main" part sets conditions for time level t+dt instead of level t,
	44	! outflow boundary conditions changed from Neumann to radiation condition,
	45	! uxrp, vynp eliminated, moisture renamed humidity
	46	!
[39]	47	! 19 2007-02-23 04:53:48Z raasch
	48	! Boundary conditions for e(nzt), pt(nzt), and q(nzt) removed because these
	49	! gridpoints are now calculated by the prognostic equation,
	50	! Dirichlet and zero gradient condition for pt established at top boundary
	51	!
[3]	52	! RCS Log replace by Id keyword, revision history cleaned up
	53	!
[1]	54	! Revision 1.15 2006/02/23 09:54:55 raasch
	55	! Surface boundary conditions in case of topography: nzb replaced by
	56	! 2d-k-index-arrays (nzb_w_inner, etc.). Conditions for u and v remain
	57	! unchanged (still using nzb) because a non-flat topography must use a
	58	! Prandtl-layer, which don't requires explicit setting of the surface values.
	59	!
	60	! Revision 1.1 1997/09/12 06:21:34 raasch
	61	! Initial revision
	62	!
	63	!
	64	! Description:
	65	! ------------
	66	! Boundary conditions for the prognostic quantities (range='main').
	67	! In case of non-cyclic lateral boundaries the conditions for velocities at
	68	! the outflow are set after the pressure solver has been called (range=
	69	! 'outflow_uvw').
	70	! One additional bottom boundary condition is applied for the TKE (=(u)*2)
	71	! in prandtl_fluxes. The cyclic lateral boundary conditions are implicitly
	72	! handled in routine exchange_horiz. Pressure boundary conditions are
	73	! explicitly set in routines pres, poisfft, poismg and sor.
	74	!------------------------------------------------------------------------------!
	75
	76	USE arrays_3d
	77	USE control_parameters
	78	USE grid_variables
	79	USE indices
	80	USE pegrid
	81
	82	IMPLICIT NONE
	83
	84	CHARACTER (LEN=*) :: range
	85
	86	INTEGER :: i, j, k
	87
[106]	88	REAL :: c_max, denom
[1]	89
[73]	90
[1]	91	IF ( range == 'main') THEN
	92	!
[667]	93	!-- Bottom boundary
	94	IF ( ibc_uv_b == 1 ) THEN
[73]	95	u_p(nzb,:,:) = u_p(nzb+1,:,:)
	96	v_p(nzb,:,:) = v_p(nzb+1,:,:)
[1]	97	ENDIF
[667]	98	DO i = nxlg, nxrg
	99	DO j = nysg, nyng
[73]	100	w_p(nzb_w_inner(j,i),j,i) = 0.0
[1]	101	ENDDO
	102	ENDDO
	103
	104	!
	105	!-- Top boundary
	106	IF ( ibc_uv_t == 0 ) THEN
[767]	107	u_p(nzt+1,:,:) = u_init(nzt+1)
	108	v_p(nzt+1,:,:) = v_init(nzt+1)
[1]	109	ELSE
[667]	110	u_p(nzt+1,:,:) = u_p(nzt,:,:)
	111	v_p(nzt+1,:,:) = v_p(nzt,:,:)
[1]	112	ENDIF
[73]	113	w_p(nzt:nzt+1,:,:) = 0.0 ! nzt is not a prognostic level (but cf. pres)
[1]	114
	115	!
[102]	116	!-- Temperature at bottom boundary.
	117	!-- In case of coupled runs (ibc_pt_b = 2) the temperature is given by
	118	!-- the sea surface temperature of the coupled ocean model.
[1]	119	IF ( ibc_pt_b == 0 ) THEN
[667]	120	DO i = nxlg, nxrg
	121	DO j = nysg, nyng
[73]	122	pt_p(nzb_s_inner(j,i),j,i) = pt(nzb_s_inner(j,i),j,i)
[1]	123	ENDDO
[73]	124	ENDDO
[102]	125	ELSEIF ( ibc_pt_b == 1 ) THEN
[667]	126	DO i = nxlg, nxrg
	127	DO j = nysg, nyng
[73]	128	pt_p(nzb_s_inner(j,i),j,i) = pt_p(nzb_s_inner(j,i)+1,j,i)
[1]	129	ENDDO
	130	ENDDO
	131	ENDIF
	132
	133	!
	134	!-- Temperature at top boundary
[19]	135	IF ( ibc_pt_t == 0 ) THEN
[667]	136	pt_p(nzt+1,:,:) = pt(nzt+1,:,:)
[19]	137	ELSEIF ( ibc_pt_t == 1 ) THEN
[667]	138	pt_p(nzt+1,:,:) = pt_p(nzt,:,:)
[19]	139	ELSEIF ( ibc_pt_t == 2 ) THEN
[667]	140	pt_p(nzt+1,:,:) = pt_p(nzt,:,:) + bc_pt_t_val * dzu(nzt+1)
[1]	141	ENDIF
	142
	143	!
	144	!-- Boundary conditions for TKE
	145	!-- Generally Neumann conditions with de/dz=0 are assumed
	146	IF ( .NOT. constant_diffusion ) THEN
[667]	147	DO i = nxlg, nxrg
	148	DO j = nysg, nyng
[73]	149	e_p(nzb_s_inner(j,i),j,i) = e_p(nzb_s_inner(j,i)+1,j,i)
[1]	150	ENDDO
	151	ENDDO
[73]	152	e_p(nzt+1,:,:) = e_p(nzt,:,:)
[1]	153	ENDIF
	154
	155	!
[95]	156	!-- Boundary conditions for salinity
	157	IF ( ocean ) THEN
	158	!
	159	!-- Bottom boundary: Neumann condition because salinity flux is always
	160	!-- given
[667]	161	DO i = nxlg, nxrg
	162	DO j = nysg, nyng
[95]	163	sa_p(nzb_s_inner(j,i),j,i) = sa_p(nzb_s_inner(j,i)+1,j,i)
	164	ENDDO
	165	ENDDO
	166
	167	!
	168	!-- Top boundary: Dirichlet or Neumann
	169	IF ( ibc_sa_t == 0 ) THEN
[667]	170	sa_p(nzt+1,:,:) = sa(nzt+1,:,:)
[95]	171	ELSEIF ( ibc_sa_t == 1 ) THEN
[667]	172	sa_p(nzt+1,:,:) = sa_p(nzt,:,:)
[95]	173	ENDIF
	174
	175	ENDIF
	176
	177	!
[1]	178	!-- Boundary conditions for total water content or scalar,
[95]	179	!-- bottom and top boundary (see also temperature)
[75]	180	IF ( humidity .OR. passive_scalar ) THEN
[1]	181	!
[75]	182	!-- Surface conditions for constant_humidity_flux
[1]	183	IF ( ibc_q_b == 0 ) THEN
[667]	184	DO i = nxlg, nxrg
	185	DO j = nysg, nyng
[73]	186	q_p(nzb_s_inner(j,i),j,i) = q(nzb_s_inner(j,i),j,i)
[1]	187	ENDDO
[73]	188	ENDDO
[1]	189	ELSE
[667]	190	DO i = nxlg, nxrg
	191	DO j = nysg, nyng
[73]	192	q_p(nzb_s_inner(j,i),j,i) = q_p(nzb_s_inner(j,i)+1,j,i)
[1]	193	ENDDO
	194	ENDDO
	195	ENDIF
	196	!
	197	!-- Top boundary
[73]	198	q_p(nzt+1,:,:) = q_p(nzt,:,:) + bc_q_t_val * dzu(nzt+1)
[667]	199
[1]	200	ENDIF
	201
	202	!
[875]	203	!-- In case of inflow at the south boundary the boundary for v is at nys
	204	!-- and in case of inflow at the left boundary the boundary for u is at nxl.
	205	!-- Since in prognostic_equations (cache optimized version) these levels are
	206	!-- handled as a prognostic level, boundary values have to be restored here.
[978]	207	!-- For the SGS-TKE, Neumann boundary conditions are used at the inflow.
[1]	208	IF ( inflow_s ) THEN
[73]	209	v_p(:,nys,:) = v_p(:,nys-1,:)
[978]	210	IF ( .NOT. constant_diffusion ) e_p(:,nys-1,:) = e_p(:,nys,:)
	211	ELSEIF ( inflow_n ) THEN
	212	IF ( .NOT. constant_diffusion ) e_p(:,nyn+1,:) = e_p(:,nyn,:)
[1]	213	ELSEIF ( inflow_l ) THEN
[73]	214	u_p(:,:,nxl) = u_p(:,:,nxl-1)
[978]	215	IF ( .NOT. constant_diffusion ) e_p(:,:,nxl-1) = e_p(:,:,nxl)
	216	ELSEIF ( inflow_r ) THEN
	217	IF ( .NOT. constant_diffusion ) e_p(:,:,nxr+1) = e_p(:,:,nxr)
[1]	218	ENDIF
	219
	220	!
	221	!-- Lateral boundary conditions for scalar quantities at the outflow
	222	IF ( outflow_s ) THEN
[73]	223	pt_p(:,nys-1,:) = pt_p(:,nys,:)
	224	IF ( .NOT. constant_diffusion ) e_p(:,nys-1,:) = e_p(:,nys,:)
[75]	225	IF ( humidity .OR. passive_scalar ) q_p(:,nys-1,:) = q_p(:,nys,:)
[1]	226	ELSEIF ( outflow_n ) THEN
[73]	227	pt_p(:,nyn+1,:) = pt_p(:,nyn,:)
	228	IF ( .NOT. constant_diffusion ) e_p(:,nyn+1,:) = e_p(:,nyn,:)
[75]	229	IF ( humidity .OR. passive_scalar ) q_p(:,nyn+1,:) = q_p(:,nyn,:)
[1]	230	ELSEIF ( outflow_l ) THEN
[73]	231	pt_p(:,:,nxl-1) = pt_p(:,:,nxl)
	232	IF ( .NOT. constant_diffusion ) e_p(:,:,nxl-1) = e_p(:,:,nxl)
[75]	233	IF ( humidity .OR. passive_scalar ) q_p(:,:,nxl-1) = q_p(:,:,nxl)
[1]	234	ELSEIF ( outflow_r ) THEN
[73]	235	pt_p(:,:,nxr+1) = pt_p(:,:,nxr)
	236	IF ( .NOT. constant_diffusion ) e_p(:,:,nxr+1) = e_p(:,:,nxr)
[75]	237	IF ( humidity .OR. passive_scalar ) q_p(:,:,nxr+1) = q_p(:,:,nxr)
[1]	238	ENDIF
	239
	240	ENDIF
	241
	242	!
[978]	243	!-- Neumann or Radiation boundary condition for the velocities at the
	244	!-- respective outflow
[106]	245	IF ( outflow_s ) THEN
[75]	246
[978]	247	IF ( bc_ns_dirneu ) THEN
	248	u(:,-1,:) = u(:,0,:)
	249	v(:,0,:) = v(:,1,:)
	250	w(:,-1,:) = w(:,0,:)
	251	ELSEIF ( bc_ns_dirrad ) THEN
[75]	252
[978]	253	c_max = dy / dt_3d
[75]	254
[978]	255	c_u_m_l = 0.0
	256	c_v_m_l = 0.0
	257	c_w_m_l = 0.0
	258
	259	c_u_m = 0.0
	260	c_v_m = 0.0
	261	c_w_m = 0.0
	262
[75]	263	!
[996]	264	!-- Calculate the phase speeds for u, v, and w, first local and then
	265	!-- average along the outflow boundary.
	266	DO k = nzb+1, nzt+1
	267	DO i = nxl, nxr
[75]	268
[106]	269	denom = u_m_s(k,0,i) - u_m_s(k,1,i)
	270
	271	IF ( denom /= 0.0 ) THEN
[996]	272	c_u(k,i) = -c_max * ( u(k,0,i) - u_m_s(k,0,i) ) / ( denom * tsc(2) )
[106]	273	IF ( c_u(k,i) < 0.0 ) THEN
	274	c_u(k,i) = 0.0
	275	ELSEIF ( c_u(k,i) > c_max ) THEN
	276	c_u(k,i) = c_max
	277	ENDIF
	278	ELSE
	279	c_u(k,i) = c_max
[75]	280	ENDIF
	281
[106]	282	denom = v_m_s(k,1,i) - v_m_s(k,2,i)
	283
	284	IF ( denom /= 0.0 ) THEN
[996]	285	c_v(k,i) = -c_max * ( v(k,1,i) - v_m_s(k,1,i) ) / ( denom * tsc(2) )
[106]	286	IF ( c_v(k,i) < 0.0 ) THEN
	287	c_v(k,i) = 0.0
	288	ELSEIF ( c_v(k,i) > c_max ) THEN
	289	c_v(k,i) = c_max
	290	ENDIF
	291	ELSE
	292	c_v(k,i) = c_max
[75]	293	ENDIF
	294
[106]	295	denom = w_m_s(k,0,i) - w_m_s(k,1,i)
[75]	296
[106]	297	IF ( denom /= 0.0 ) THEN
[996]	298	c_w(k,i) = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) / ( denom * tsc(2) )
[106]	299	IF ( c_w(k,i) < 0.0 ) THEN
	300	c_w(k,i) = 0.0
	301	ELSEIF ( c_w(k,i) > c_max ) THEN
	302	c_w(k,i) = c_max
	303	ENDIF
	304	ELSE
	305	c_w(k,i) = c_max
[75]	306	ENDIF
[106]	307
[978]	308	c_u_m_l(k) = c_u_m_l(k) + c_u(k,i)
	309	c_v_m_l(k) = c_v_m_l(k) + c_v(k,i)
	310	c_w_m_l(k) = c_w_m_l(k) + c_w(k,i)
[106]	311
[978]	312	ENDDO
	313	ENDDO
[75]	314
[978]	315	#if defined( __parallel )
	316	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
	317	CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
	318	MPI_SUM, comm1dx, ierr )
	319	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
	320	CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
	321	MPI_SUM, comm1dx, ierr )
	322	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
	323	CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
	324	MPI_SUM, comm1dx, ierr )
	325	#else
	326	c_u_m = c_u_m_l
	327	c_v_m = c_v_m_l
	328	c_w_m = c_w_m_l
	329	#endif
	330
	331	c_u_m = c_u_m / (nx+1)
	332	c_v_m = c_v_m / (nx+1)
	333	c_w_m = c_w_m / (nx+1)
	334
[75]	335	!
[978]	336	!-- Save old timelevels for the next timestep
	337	IF ( intermediate_timestep_count == 1 ) THEN
	338	u_m_s(:,:,:) = u(:,0:1,:)
	339	v_m_s(:,:,:) = v(:,1:2,:)
	340	w_m_s(:,:,:) = w(:,0:1,:)
	341	ENDIF
	342
	343	!
	344	!-- Calculate the new velocities
[996]	345	DO k = nzb+1, nzt+1
	346	DO i = nxlg, nxrg
[978]	347	u_p(k,-1,i) = u(k,-1,i) - dt_3d * tsc(2) * c_u_m(k) * &
[75]	348	( u(k,-1,i) - u(k,0,i) ) * ddy
	349
[978]	350	v_p(k,0,i) = v(k,0,i) - dt_3d * tsc(2) * c_v_m(k) * &
[106]	351	( v(k,0,i) - v(k,1,i) ) * ddy
[75]	352
[978]	353	w_p(k,-1,i) = w(k,-1,i) - dt_3d * tsc(2) * c_w_m(k) * &
[75]	354	( w(k,-1,i) - w(k,0,i) ) * ddy
[978]	355	ENDDO
[75]	356	ENDDO
	357
	358	!
[978]	359	!-- Bottom boundary at the outflow
	360	IF ( ibc_uv_b == 0 ) THEN
	361	u_p(nzb,-1,:) = 0.0
	362	v_p(nzb,0,:) = 0.0
	363	ELSE
	364	u_p(nzb,-1,:) = u_p(nzb+1,-1,:)
	365	v_p(nzb,0,:) = v_p(nzb+1,0,:)
	366	ENDIF
	367	w_p(nzb,-1,:) = 0.0
[73]	368
[75]	369	!
[978]	370	!-- Top boundary at the outflow
	371	IF ( ibc_uv_t == 0 ) THEN
	372	u_p(nzt+1,-1,:) = u_init(nzt+1)
	373	v_p(nzt+1,0,:) = v_init(nzt+1)
	374	ELSE
	375	u_p(nzt+1,-1,:) = u(nzt,-1,:)
	376	v_p(nzt+1,0,:) = v(nzt,0,:)
	377	ENDIF
	378	w_p(nzt:nzt+1,-1,:) = 0.0
	379
[75]	380	ENDIF
[73]	381
[75]	382	ENDIF
[73]	383
[106]	384	IF ( outflow_n ) THEN
[73]	385
[978]	386	IF ( bc_ns_neudir ) THEN
	387	u(:,ny+1,:) = u(:,ny,:)
	388	v(:,ny+1,:) = v(:,ny,:)
	389	w(:,ny+1,:) = w(:,ny,:)
	390	ELSEIF ( bc_ns_dirrad ) THEN
[75]	391
[978]	392	c_max = dy / dt_3d
[75]	393
[978]	394	c_u_m_l = 0.0
	395	c_v_m_l = 0.0
	396	c_w_m_l = 0.0
	397
	398	c_u_m = 0.0
	399	c_v_m = 0.0
	400	c_w_m = 0.0
	401
[1]	402	!
[996]	403	!-- Calculate the phase speeds for u, v, and w, first local and then
	404	!-- average along the outflow boundary.
	405	DO k = nzb+1, nzt+1
	406	DO i = nxl, nxr
[73]	407
[106]	408	denom = u_m_n(k,ny,i) - u_m_n(k,ny-1,i)
	409
	410	IF ( denom /= 0.0 ) THEN
[996]	411	c_u(k,i) = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) / ( denom * tsc(2) )
[106]	412	IF ( c_u(k,i) < 0.0 ) THEN
	413	c_u(k,i) = 0.0
	414	ELSEIF ( c_u(k,i) > c_max ) THEN
	415	c_u(k,i) = c_max
	416	ENDIF
	417	ELSE
	418	c_u(k,i) = c_max
[73]	419	ENDIF
	420
[106]	421	denom = v_m_n(k,ny,i) - v_m_n(k,ny-1,i)
[73]	422
[106]	423	IF ( denom /= 0.0 ) THEN
[996]	424	c_v(k,i) = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / ( denom * tsc(2) )
[106]	425	IF ( c_v(k,i) < 0.0 ) THEN
	426	c_v(k,i) = 0.0
	427	ELSEIF ( c_v(k,i) > c_max ) THEN
	428	c_v(k,i) = c_max
	429	ENDIF
	430	ELSE
	431	c_v(k,i) = c_max
[73]	432	ENDIF
	433
[106]	434	denom = w_m_n(k,ny,i) - w_m_n(k,ny-1,i)
[73]	435
[106]	436	IF ( denom /= 0.0 ) THEN
[996]	437	c_w(k,i) = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / ( denom * tsc(2) )
[106]	438	IF ( c_w(k,i) < 0.0 ) THEN
	439	c_w(k,i) = 0.0
	440	ELSEIF ( c_w(k,i) > c_max ) THEN
	441	c_w(k,i) = c_max
	442	ENDIF
	443	ELSE
	444	c_w(k,i) = c_max
[73]	445	ENDIF
[106]	446
[978]	447	c_u_m_l(k) = c_u_m_l(k) + c_u(k,i)
	448	c_v_m_l(k) = c_v_m_l(k) + c_v(k,i)
	449	c_w_m_l(k) = c_w_m_l(k) + c_w(k,i)
[106]	450
[978]	451	ENDDO
	452	ENDDO
[73]	453
[978]	454	#if defined( __parallel )
	455	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
	456	CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
	457	MPI_SUM, comm1dx, ierr )
	458	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
	459	CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
	460	MPI_SUM, comm1dx, ierr )
	461	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
	462	CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
	463	MPI_SUM, comm1dx, ierr )
	464	#else
	465	c_u_m = c_u_m_l
	466	c_v_m = c_v_m_l
	467	c_w_m = c_w_m_l
	468	#endif
	469
	470	c_u_m = c_u_m / (nx+1)
	471	c_v_m = c_v_m / (nx+1)
	472	c_w_m = c_w_m / (nx+1)
	473
[73]	474	!
[978]	475	!-- Save old timelevels for the next timestep
	476	IF ( intermediate_timestep_count == 1 ) THEN
	477	u_m_n(:,:,:) = u(:,ny-1:ny,:)
	478	v_m_n(:,:,:) = v(:,ny-1:ny,:)
	479	w_m_n(:,:,:) = w(:,ny-1:ny,:)
	480	ENDIF
[73]	481
[978]	482	!
	483	!-- Calculate the new velocities
[996]	484	DO k = nzb+1, nzt+1
	485	DO i = nxlg, nxrg
[978]	486	u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * tsc(2) * c_u_m(k) * &
	487	( u(k,ny+1,i) - u(k,ny,i) ) * ddy
[73]	488
[978]	489	v_p(k,ny+1,i) = v(k,ny+1,i) - dt_3d * tsc(2) * c_v_m(k) * &
	490	( v(k,ny+1,i) - v(k,ny,i) ) * ddy
[73]	491
[978]	492	w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * tsc(2) * c_w_m(k) * &
	493	( w(k,ny+1,i) - w(k,ny,i) ) * ddy
	494	ENDDO
[1]	495	ENDDO
	496
	497	!
[978]	498	!-- Bottom boundary at the outflow
	499	IF ( ibc_uv_b == 0 ) THEN
	500	u_p(nzb,ny+1,:) = 0.0
	501	v_p(nzb,ny+1,:) = 0.0
	502	ELSE
	503	u_p(nzb,ny+1,:) = u_p(nzb+1,ny+1,:)
	504	v_p(nzb,ny+1,:) = v_p(nzb+1,ny+1,:)
	505	ENDIF
	506	w_p(nzb,ny+1,:) = 0.0
[73]	507
	508	!
[978]	509	!-- Top boundary at the outflow
	510	IF ( ibc_uv_t == 0 ) THEN
	511	u_p(nzt+1,ny+1,:) = u_init(nzt+1)
	512	v_p(nzt+1,ny+1,:) = v_init(nzt+1)
	513	ELSE
	514	u_p(nzt+1,ny+1,:) = u_p(nzt,nyn+1,:)
	515	v_p(nzt+1,ny+1,:) = v_p(nzt,nyn+1,:)
	516	ENDIF
	517	w_p(nzt:nzt+1,ny+1,:) = 0.0
	518
[1]	519	ENDIF
	520
[75]	521	ENDIF
	522
[106]	523	IF ( outflow_l ) THEN
[75]	524
[978]	525	IF ( bc_lr_neudir ) THEN
	526	u(:,:,-1) = u(:,:,0)
	527	v(:,:,0) = v(:,:,1)
	528	w(:,:,-1) = w(:,:,0)
	529	ELSEIF ( bc_ns_dirrad ) THEN
[75]	530
[978]	531	c_max = dx / dt_3d
[75]	532
[978]	533	c_u_m_l = 0.0
	534	c_v_m_l = 0.0
	535	c_w_m_l = 0.0
	536
	537	c_u_m = 0.0
	538	c_v_m = 0.0
	539	c_w_m = 0.0
	540
[1]	541	!
[996]	542	!-- Calculate the phase speeds for u, v, and w, first local and then
	543	!-- average along the outflow boundary.
	544	DO k = nzb+1, nzt+1
	545	DO j = nys, nyn
[75]	546
[106]	547	denom = u_m_l(k,j,1) - u_m_l(k,j,2)
	548
	549	IF ( denom /= 0.0 ) THEN
[996]	550	c_u(k,j) = -c_max * ( u(k,j,1) - u_m_l(k,j,1) ) / ( denom * tsc(2) )
[107]	551	IF ( c_u(k,j) < 0.0 ) THEN
[106]	552	c_u(k,j) = 0.0
[107]	553	ELSEIF ( c_u(k,j) > c_max ) THEN
	554	c_u(k,j) = c_max
[106]	555	ENDIF
	556	ELSE
[107]	557	c_u(k,j) = c_max
[75]	558	ENDIF
	559
[106]	560	denom = v_m_l(k,j,0) - v_m_l(k,j,1)
[75]	561
[106]	562	IF ( denom /= 0.0 ) THEN
[996]	563	c_v(k,j) = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) / ( denom * tsc(2) )
[106]	564	IF ( c_v(k,j) < 0.0 ) THEN
	565	c_v(k,j) = 0.0
	566	ELSEIF ( c_v(k,j) > c_max ) THEN
	567	c_v(k,j) = c_max
	568	ENDIF
	569	ELSE
	570	c_v(k,j) = c_max
[75]	571	ENDIF
	572
[106]	573	denom = w_m_l(k,j,0) - w_m_l(k,j,1)
[75]	574
[106]	575	IF ( denom /= 0.0 ) THEN
[996]	576	c_w(k,j) = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) / ( denom * tsc(2) )
[106]	577	IF ( c_w(k,j) < 0.0 ) THEN
	578	c_w(k,j) = 0.0
	579	ELSEIF ( c_w(k,j) > c_max ) THEN
	580	c_w(k,j) = c_max
	581	ENDIF
	582	ELSE
	583	c_w(k,j) = c_max
[75]	584	ENDIF
[106]	585
[978]	586	c_u_m_l(k) = c_u_m_l(k) + c_u(k,j)
	587	c_v_m_l(k) = c_v_m_l(k) + c_v(k,j)
	588	c_w_m_l(k) = c_w_m_l(k) + c_w(k,j)
[106]	589
[978]	590	ENDDO
	591	ENDDO
[75]	592
[978]	593	#if defined( __parallel )
	594	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
	595	CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
	596	MPI_SUM, comm1dy, ierr )
	597	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
	598	CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
	599	MPI_SUM, comm1dy, ierr )
	600	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
	601	CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
	602	MPI_SUM, comm1dy, ierr )
	603	#else
	604	c_u_m = c_u_m_l
	605	c_v_m = c_v_m_l
	606	c_w_m = c_w_m_l
	607	#endif
	608
	609	c_u_m = c_u_m / (ny+1)
	610	c_v_m = c_v_m / (ny+1)
	611	c_w_m = c_w_m / (ny+1)
	612
[73]	613	!
[978]	614	!-- Save old timelevels for the next timestep
	615	IF ( intermediate_timestep_count == 1 ) THEN
	616	u_m_l(:,:,:) = u(:,:,1:2)
	617	v_m_l(:,:,:) = v(:,:,0:1)
	618	w_m_l(:,:,:) = w(:,:,0:1)
	619	ENDIF
	620
	621	!
	622	!-- Calculate the new velocities
[996]	623	DO k = nzb+1, nzt+1
	624	DO i = nxlg, nxrg
[978]	625	u_p(k,j,0) = u(k,j,0) - dt_3d * tsc(2) * c_u_m(k) * &
[106]	626	( u(k,j,0) - u(k,j,1) ) * ddx
[75]	627
[978]	628	v_p(k,j,-1) = v(k,j,-1) - dt_3d * tsc(2) * c_v_m(k) * &
[75]	629	( v(k,j,-1) - v(k,j,0) ) * ddx
	630
[978]	631	w_p(k,j,-1) = w(k,j,-1) - dt_3d * tsc(2) * c_w_m(k) * &
[75]	632	( w(k,j,-1) - w(k,j,0) ) * ddx
[978]	633	ENDDO
[75]	634	ENDDO
	635
	636	!
[978]	637	!-- Bottom boundary at the outflow
	638	IF ( ibc_uv_b == 0 ) THEN
	639	u_p(nzb,:,0) = 0.0
	640	v_p(nzb,:,-1) = 0.0
	641	ELSE
	642	u_p(nzb,:,0) = u_p(nzb+1,:,0)
	643	v_p(nzb,:,-1) = v_p(nzb+1,:,-1)
	644	ENDIF
	645	w_p(nzb,:,-1) = 0.0
[1]	646
[75]	647	!
[978]	648	!-- Top boundary at the outflow
	649	IF ( ibc_uv_t == 0 ) THEN
	650	u_p(nzt+1,:,-1) = u_init(nzt+1)
	651	v_p(nzt+1,:,-1) = v_init(nzt+1)
	652	ELSE
	653	u_p(nzt+1,:,-1) = u_p(nzt,:,-1)
	654	v_p(nzt+1,:,-1) = v_p(nzt,:,-1)
	655	ENDIF
	656	w_p(nzt:nzt+1,:,-1) = 0.0
	657
[75]	658	ENDIF
[73]	659
[75]	660	ENDIF
[73]	661
[106]	662	IF ( outflow_r ) THEN
[73]	663
[978]	664	IF ( bc_lr_dirneu ) THEN
	665	u(:,:,nx+1) = u(:,:,nx)
	666	v(:,:,nx+1) = v(:,:,nx)
	667	w(:,:,nx+1) = w(:,:,nx)
	668	ELSEIF ( bc_ns_dirrad ) THEN
[75]	669
[978]	670	c_max = dx / dt_3d
[75]	671
[978]	672	c_u_m_l = 0.0
	673	c_v_m_l = 0.0
	674	c_w_m_l = 0.0
	675
	676	c_u_m = 0.0
	677	c_v_m = 0.0
	678	c_w_m = 0.0
	679
[1]	680	!
[996]	681	!-- Calculate the phase speeds for u, v, and w, first local and then
	682	!-- average along the outflow boundary.
	683	DO k = nzb+1, nzt+1
	684	DO j = nys, nyn
[73]	685
[106]	686	denom = u_m_r(k,j,nx) - u_m_r(k,j,nx-1)
	687
	688	IF ( denom /= 0.0 ) THEN
[996]	689	c_u(k,j) = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) ) / ( denom * tsc(2) )
[106]	690	IF ( c_u(k,j) < 0.0 ) THEN
	691	c_u(k,j) = 0.0
	692	ELSEIF ( c_u(k,j) > c_max ) THEN
	693	c_u(k,j) = c_max
	694	ENDIF
	695	ELSE
	696	c_u(k,j) = c_max
[73]	697	ENDIF
	698
[106]	699	denom = v_m_r(k,j,nx) - v_m_r(k,j,nx-1)
[73]	700
[106]	701	IF ( denom /= 0.0 ) THEN
[996]	702	c_v(k,j) = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) / ( denom * tsc(2) )
[106]	703	IF ( c_v(k,j) < 0.0 ) THEN
	704	c_v(k,j) = 0.0
	705	ELSEIF ( c_v(k,j) > c_max ) THEN
	706	c_v(k,j) = c_max
	707	ENDIF
	708	ELSE
	709	c_v(k,j) = c_max
[73]	710	ENDIF
	711
[106]	712	denom = w_m_r(k,j,nx) - w_m_r(k,j,nx-1)
[73]	713
[106]	714	IF ( denom /= 0.0 ) THEN
[996]	715	c_w(k,j) = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) / ( denom * tsc(2) )
[106]	716	IF ( c_w(k,j) < 0.0 ) THEN
	717	c_w(k,j) = 0.0
	718	ELSEIF ( c_w(k,j) > c_max ) THEN
	719	c_w(k,j) = c_max
	720	ENDIF
	721	ELSE
	722	c_w(k,j) = c_max
[73]	723	ENDIF
[106]	724
[978]	725	c_u_m_l(k) = c_u_m_l(k) + c_u(k,j)
	726	c_v_m_l(k) = c_v_m_l(k) + c_v(k,j)
	727	c_w_m_l(k) = c_w_m_l(k) + c_w(k,j)
[106]	728
[978]	729	ENDDO
	730	ENDDO
[73]	731
[978]	732	#if defined( __parallel )
	733	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
	734	CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
	735	MPI_SUM, comm1dy, ierr )
	736	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
	737	CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
	738	MPI_SUM, comm1dy, ierr )
	739	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
	740	CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
	741	MPI_SUM, comm1dy, ierr )
	742	#else
	743	c_u_m = c_u_m_l
	744	c_v_m = c_v_m_l
	745	c_w_m = c_w_m_l
	746	#endif
	747
	748	c_u_m = c_u_m / (ny+1)
	749	c_v_m = c_v_m / (ny+1)
	750	c_w_m = c_w_m / (ny+1)
	751
[73]	752	!
[978]	753	!-- Save old timelevels for the next timestep
	754	IF ( intermediate_timestep_count == 1 ) THEN
	755	u_m_r(:,:,:) = u(:,:,nx-1:nx)
	756	v_m_r(:,:,:) = v(:,:,nx-1:nx)
	757	w_m_r(:,:,:) = w(:,:,nx-1:nx)
	758	ENDIF
[73]	759
[978]	760	!
	761	!-- Calculate the new velocities
[996]	762	DO k = nzb+1, nzt+1
	763	DO i = nxlg, nxrg
[978]	764	u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * tsc(2) * c_u_m(k) * &
	765	( u(k,j,nx+1) - u(k,j,nx) ) * ddx
[73]	766
[978]	767	v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * tsc(2) * c_v_m(k) * &
	768	( v(k,j,nx+1) - v(k,j,nx) ) * ddx
[73]	769
[978]	770	w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * tsc(2) * c_w_m(k) * &
	771	( w(k,j,nx+1) - w(k,j,nx) ) * ddx
	772	ENDDO
[73]	773	ENDDO
	774
	775	!
[978]	776	!-- Bottom boundary at the outflow
	777	IF ( ibc_uv_b == 0 ) THEN
	778	u_p(nzb,:,nx+1) = 0.0
	779	v_p(nzb,:,nx+1) = 0.0
	780	ELSE
	781	u_p(nzb,:,nx+1) = u_p(nzb+1,:,nx+1)
	782	v_p(nzb,:,nx+1) = v_p(nzb+1,:,nx+1)
	783	ENDIF
	784	w_p(nzb,:,nx+1) = 0.0
[73]	785
	786	!
[978]	787	!-- Top boundary at the outflow
	788	IF ( ibc_uv_t == 0 ) THEN
	789	u_p(nzt+1,:,nx+1) = u_init(nzt+1)
	790	v_p(nzt+1,:,nx+1) = v_init(nzt+1)
	791	ELSE
	792	u_p(nzt+1,:,nx+1) = u_p(nzt,:,nx+1)
	793	v_p(nzt+1,:,nx+1) = v_p(nzt,:,nx+1)
	794	ENDIF
	795	w(nzt:nzt+1,:,nx+1) = 0.0
	796
[1]	797	ENDIF
	798
	799	ENDIF
	800
	801	END SUBROUTINE boundary_conds

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |