Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |