Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |