Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |