Home

Context Navigation

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

Last change on this file since 1339 was 1321, checked in by raasch, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 30.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |