Home

Context Navigation

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

Last change on this file since 1093 was 1054, checked in by hoffmann, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 29.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |