Home

Context Navigation

boundary_conds.f90 @ 91

Last change on this file since 91 was 77, checked in by raasch, 18 years ago

New:
---

particle reflection from vertical walls implemented, particle SGS model adjusted to walls

Wall functions for vertical walls now include diabatic conditions. New subroutines wall_fluxes, wall_fluxes_e. New 4D-array rif_wall.

new d3par-parameter netcdf_64bit_3d to switch on 64bit offset only for 3D files

new d3par-parameter dt_max to define the maximum value for the allowed timestep

new inipar-parameter loop_optimization to control the loop optimization method

new inipar-parameter pt_refrence. If given, this value is used as the reference that in buoyancy terms (otherwise, the instantaneous horizontally averaged temperature is used).

new user interface user_advec_particles

new initializing action "by_user" calls user_init_3d_model and allows the initial setting of all 3d arrays

topography height informations are stored on arrays zu_s_inner and zw_w_inner and output to the 2d/3d NetCDF files

samples added to the user interface which show how to add user-define time series quantities.

calculation/output of precipitation amount, precipitation rate and z0 (by setting "pra*", "prr*", "z0*" with data_output). The time interval on which the precipitation amount is defined is set by new d3par-parameter precipitation_amount_interval

unit 9 opened for debug output (file DEBUG_<pe#>)

Makefile, advec_particles, average_3d_data, buoyancy, calc_precipitation, check_open, check_parameters, data_output_2d, diffusion_e, diffusion_u, diffusion_v, diffusion_w, diffusivities, header, impact_of_latent_heat, init_particles, init_3d_model, modules, netcdf, parin, production_e, read_var_list, read_3d_binary, sum_up_3d_data, user_interface, write_var_list, write_3d_binary

New: wall_fluxes

Changed:

General revision of non-cyclic horizontal boundary conditions: radiation boundary conditions are now used instead of Neumann conditions at the outflow (calculation needs velocity values for t-dt, which are stored on new arrays u_m_l, u_m_r, etc.), calculation of mean outflow is not needed any more, volume flow control is added for the outflow boundary (currently only for the north boundary!!), additional gridpoints along x and y (uxrp, vynp) are not needed any more, routine "boundary_conds" now operates on timelevel t+dt and is not split in two parts (main, uvw_outflow) any more, Neumann boundary conditions at inflow/outflow in case of non-cyclic boundary conditions for all 2d-arrays that are handled by exchange_horiz_2d

The FFT-method for solving the Poisson-equation is now working with Neumann boundary conditions both at the bottom and the top. This requires adjustments of the tridiagonal coefficients and subtracting the horizontally averaged mean from the vertical velocity field.

+age_m in particle_type

Particles-package is now part of the default code ("-p particles" is not needed any more).

Move call of user_actions( 'after_integration' ) below increment of times
and counters. user_actions is now called for each statistic region and has as an argument the number of the respective region (sr)

d3par-parameter data_output_ts removed. Timeseries output for "profil" removed. Timeseries are now switched on by dt_dots. Timeseries data is collected in flow_statistics.

Initial velocities at nzb+1 are regarded for volume flow control in case they have been set zero before (to avoid small timesteps); see new internal parameters u/v_nzb_p1_for_vfc.

q is not allowed to become negative (prognostic_equations).

poisfft_init is only called if fft-solver is switched on (init_pegrid).

d3par-parameter moisture renamed to humidity.

Subversion global revision number is read from mrun and added to the run description header and to the run control (_rc) file.

vtk directives removed from main program.

The uitility routine interpret_config reads PALM environment variables from NAMELIST instead using the system call GETENV.

advec_u_pw, advec_u_up, advec_v_pw, advec_v_up, asselin_filter, check_parameters, coriolis, data_output_dvrp, data_output_ptseries, data_output_ts, data_output_2d, data_output_3d, diffusion_u, diffusion_v, exchange_horiz, exchange_horiz_2d, flow_statistics, header, init_grid, init_particles, init_pegrid, init_rankine, init_pt_anomaly, init_1d_model, init_3d_model, modules, palm, package_parin, parin, poisfft, poismg, prandtl_fluxes, pres, production_e, prognostic_equations, read_var_list, read_3d_binary, sor, swap_timelevel, time_integration, write_var_list, write_3d_binary

Errors:

Bugfix: preset of tendencies te_em, te_um, te_vm in init_1d_model

Bugfix in sample for reading user defined data from restart file (user_init)

Bugfix in setting diffusivities for cases with the outflow damping layer extending over more than one subdomain (init_3d_model)

Check for possible negative humidities in the initial humidity profile.

in Makefile, default suffixes removed from the suffix list to avoid calling of m2c in
# case of .mod files

Makefile
check_parameters, init_1d_model, init_3d_model, user_interface

Property svn:keywords set to Id

File size: 17.1 KB

Rev	Line
[1]	1	SUBROUTINE boundary_conds( range )
	2
	3	!------------------------------------------------------------------------------!
	4	! Actual revisions:
	5	! -----------------
[77]	6	!
[1]	7	!
	8	! Former revisions:
	9	! -----------------
[3]	10	! $Id: boundary_conds.f90 77 2007-03-29 04:26:56Z raasch $
[39]	11	!
[77]	12	! 75 2007-03-22 09:54:05Z raasch
	13	! The "main" part sets conditions for time level t+dt instead of level t,
	14	! outflow boundary conditions changed from Neumann to radiation condition,
	15	! uxrp, vynp eliminated, moisture renamed humidity
	16	!
[39]	17	! 19 2007-02-23 04:53:48Z raasch
	18	! Boundary conditions for e(nzt), pt(nzt), and q(nzt) removed because these
	19	! gridpoints are now calculated by the prognostic equation,
	20	! Dirichlet and zero gradient condition for pt established at top boundary
	21	!
[3]	22	! RCS Log replace by Id keyword, revision history cleaned up
	23	!
[1]	24	! Revision 1.15 2006/02/23 09:54:55 raasch
	25	! Surface boundary conditions in case of topography: nzb replaced by
	26	! 2d-k-index-arrays (nzb_w_inner, etc.). Conditions for u and v remain
	27	! unchanged (still using nzb) because a non-flat topography must use a
	28	! Prandtl-layer, which don't requires explicit setting of the surface values.
	29	!
	30	! Revision 1.1 1997/09/12 06:21:34 raasch
	31	! Initial revision
	32	!
	33	!
	34	! Description:
	35	! ------------
	36	! Boundary conditions for the prognostic quantities (range='main').
	37	! In case of non-cyclic lateral boundaries the conditions for velocities at
	38	! the outflow are set after the pressure solver has been called (range=
	39	! 'outflow_uvw').
	40	! One additional bottom boundary condition is applied for the TKE (=(u)*2)
	41	! in prandtl_fluxes. The cyclic lateral boundary conditions are implicitly
	42	! handled in routine exchange_horiz. Pressure boundary conditions are
	43	! explicitly set in routines pres, poisfft, poismg and sor.
	44	!------------------------------------------------------------------------------!
	45
	46	USE arrays_3d
	47	USE control_parameters
	48	USE grid_variables
	49	USE indices
	50	USE pegrid
	51
	52	IMPLICIT NONE
	53
	54	CHARACTER (LEN=*) :: range
	55
	56	INTEGER :: i, j, k
	57
[73]	58	REAL :: c_max, c_u, c_v, c_w, denom
[1]	59
[73]	60
[1]	61	IF ( range == 'main') THEN
	62	!
	63	!-- Bottom boundary
	64	IF ( ibc_uv_b == 0 ) THEN
[73]	65	!
	66	!-- Satisfying the Dirichlet condition with an extra layer below the
	67	!-- surface where the u and v component change their sign
	68	u_p(nzb,:,:) = -u_p(nzb+1,:,:)
	69	v_p(nzb,:,:) = -v_p(nzb+1,:,:)
	70	ELSE
	71	u_p(nzb,:,:) = u_p(nzb+1,:,:)
	72	v_p(nzb,:,:) = v_p(nzb+1,:,:)
[1]	73	ENDIF
	74	DO i = nxl-1, nxr+1
	75	DO j = nys-1, nyn+1
[73]	76	w_p(nzb_w_inner(j,i),j,i) = 0.0
[1]	77	ENDDO
	78	ENDDO
	79
	80	!
	81	!-- Top boundary
	82	IF ( ibc_uv_t == 0 ) THEN
[73]	83	u_p(nzt+1,:,:) = ug(nzt+1)
	84	v_p(nzt+1,:,:) = vg(nzt+1)
[1]	85	ELSE
[73]	86	u_p(nzt+1,:,:) = u_p(nzt,:,:)
	87	v_p(nzt+1,:,:) = v_p(nzt,:,:)
[1]	88	ENDIF
[73]	89	w_p(nzt:nzt+1,:,:) = 0.0 ! nzt is not a prognostic level (but cf. pres)
[1]	90
	91	!
	92	!-- Temperature at bottom boundary
	93	IF ( ibc_pt_b == 0 ) THEN
[73]	94	DO i = nxl-1, nxr+1
	95	DO j = nys-1, nyn+1
	96	pt_p(nzb_s_inner(j,i),j,i) = pt(nzb_s_inner(j,i),j,i)
[1]	97	ENDDO
[73]	98	ENDDO
[1]	99	ELSE
	100	DO i = nxl-1, nxr+1
	101	DO j = nys-1, nyn+1
[73]	102	pt_p(nzb_s_inner(j,i),j,i) = pt_p(nzb_s_inner(j,i)+1,j,i)
[1]	103	ENDDO
	104	ENDDO
	105	ENDIF
	106
	107	!
	108	!-- Temperature at top boundary
[19]	109	IF ( ibc_pt_t == 0 ) THEN
[73]	110	pt_p(nzt+1,:,:) = pt(nzt+1,:,:)
[19]	111	ELSEIF ( ibc_pt_t == 1 ) THEN
[73]	112	pt_p(nzt+1,:,:) = pt_p(nzt,:,:)
[19]	113	ELSEIF ( ibc_pt_t == 2 ) THEN
[73]	114	pt_p(nzt+1,:,:) = pt_p(nzt,:,:) + bc_pt_t_val * dzu(nzt+1)
[1]	115	ENDIF
	116
	117	!
	118	!-- Boundary conditions for TKE
	119	!-- Generally Neumann conditions with de/dz=0 are assumed
	120	IF ( .NOT. constant_diffusion ) THEN
	121	DO i = nxl-1, nxr+1
	122	DO j = nys-1, nyn+1
[73]	123	e_p(nzb_s_inner(j,i),j,i) = e_p(nzb_s_inner(j,i)+1,j,i)
[1]	124	ENDDO
	125	ENDDO
[73]	126	e_p(nzt+1,:,:) = e_p(nzt,:,:)
[1]	127	ENDIF
	128
	129	!
	130	!-- Boundary conditions for total water content or scalar,
	131	!-- bottom and surface boundary (see also temperature)
[75]	132	IF ( humidity .OR. passive_scalar ) THEN
[1]	133	!
[75]	134	!-- Surface conditions for constant_humidity_flux
[1]	135	IF ( ibc_q_b == 0 ) THEN
[73]	136	DO i = nxl-1, nxr+1
	137	DO j = nys-1, nyn+1
	138	q_p(nzb_s_inner(j,i),j,i) = q(nzb_s_inner(j,i),j,i)
[1]	139	ENDDO
[73]	140	ENDDO
[1]	141	ELSE
	142	DO i = nxl-1, nxr+1
	143	DO j = nys-1, nyn+1
[73]	144	q_p(nzb_s_inner(j,i),j,i) = q_p(nzb_s_inner(j,i)+1,j,i)
[1]	145	ENDDO
	146	ENDDO
	147	ENDIF
	148	!
	149	!-- Top boundary
[73]	150	q_p(nzt+1,:,:) = q_p(nzt,:,:) + bc_q_t_val * dzu(nzt+1)
[1]	151	ENDIF
	152
	153	!
	154	!-- Lateral boundary conditions at the inflow. Quasi Neumann conditions
	155	!-- are needed for the wall normal velocity in order to ensure zero
	156	!-- divergence. Dirichlet conditions are used for all other quantities.
	157	IF ( inflow_s ) THEN
[73]	158	v_p(:,nys,:) = v_p(:,nys-1,:)
[1]	159	ELSEIF ( inflow_n ) THEN
[75]	160	v_p(:,nyn,:) = v_p(:,nyn+1,:)
[1]	161	ELSEIF ( inflow_l ) THEN
[73]	162	u_p(:,:,nxl) = u_p(:,:,nxl-1)
[1]	163	ELSEIF ( inflow_r ) THEN
[75]	164	u_p(:,:,nxr) = u_p(:,:,nxr+1)
[1]	165	ENDIF
	166
	167	!
	168	!-- Lateral boundary conditions for scalar quantities at the outflow
	169	IF ( outflow_s ) THEN
[73]	170	pt_p(:,nys-1,:) = pt_p(:,nys,:)
	171	IF ( .NOT. constant_diffusion ) e_p(:,nys-1,:) = e_p(:,nys,:)
[75]	172	IF ( humidity .OR. passive_scalar ) q_p(:,nys-1,:) = q_p(:,nys,:)
[1]	173	ELSEIF ( outflow_n ) THEN
[73]	174	pt_p(:,nyn+1,:) = pt_p(:,nyn,:)
	175	IF ( .NOT. constant_diffusion ) e_p(:,nyn+1,:) = e_p(:,nyn,:)
[75]	176	IF ( humidity .OR. passive_scalar ) q_p(:,nyn+1,:) = q_p(:,nyn,:)
[1]	177	ELSEIF ( outflow_l ) THEN
[73]	178	pt_p(:,:,nxl-1) = pt_p(:,:,nxl)
	179	IF ( .NOT. constant_diffusion ) e_p(:,:,nxl-1) = e_p(:,:,nxl)
[75]	180	IF ( humidity .OR. passive_scalar ) q_p(:,:,nxl-1) = q_p(:,:,nxl)
[1]	181	ELSEIF ( outflow_r ) THEN
[73]	182	pt_p(:,:,nxr+1) = pt_p(:,:,nxr)
	183	IF ( .NOT. constant_diffusion ) e_p(:,:,nxr+1) = e_p(:,:,nxr)
[75]	184	IF ( humidity .OR. passive_scalar ) q_p(:,:,nxr+1) = q_p(:,:,nxr)
[1]	185	ENDIF
	186
	187	ENDIF
	188
	189	!
[75]	190	!-- Radiation boundary condition for the velocities at the respective outflow
	191	IF ( outflow_s .AND. &
	192	intermediate_timestep_count == intermediate_timestep_count_max ) &
	193	THEN
	194
	195	c_max = dy / dt_3d
	196
	197	DO i = nxl-1, nxr+1
	198	DO k = nzb+1, nzt+1
	199
	200	!
	201	!-- First calculate the phase speeds for u,v, and w
	202	denom = u_m_s(k,0,i) - u_m_s(k,1,i)
	203
	204	IF ( denom /= 0.0 ) THEN
	205	c_u = -c_max * ( u(k,0,i) - u_m_s(k,0,i) ) / denom
	206	IF ( c_u < 0.0 ) THEN
	207	c_u = 0.0
	208	ELSEIF ( c_u > c_max ) THEN
	209	c_u = c_max
	210	ENDIF
	211	ELSE
	212	c_u = c_max
	213	ENDIF
	214	denom = v_m_s(k,0,i) - v_m_s(k,1,i)
	215
	216	IF ( denom /= 0.0 ) THEN
	217	c_v = -c_max * ( v(k,0,i) - v_m_s(k,0,i) ) / denom
	218	IF ( c_v < 0.0 ) THEN
	219	c_v = 0.0
	220	ELSEIF ( c_v > c_max ) THEN
	221	c_v = c_max
	222	ENDIF
	223	ELSE
	224	c_v = c_max
	225	ENDIF
	226
	227	denom = w_m_s(k,0,i) - w_m_s(k,1,i)
	228
	229	IF ( denom /= 0.0 ) THEN
	230	c_w = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) / denom
	231	IF ( c_w < 0.0 ) THEN
	232	c_w = 0.0
	233	ELSEIF ( c_w > c_max ) THEN
	234	c_w = c_max
	235	ENDIF
	236	ELSE
	237	c_w = c_max
	238	ENDIF
	239
	240	!
	241	!-- Calculate the new velocities
	242	u_p(k,-1,i) = u(k,-1,i) + dt_3d * c_u * &
	243	( u(k,-1,i) - u(k,0,i) ) * ddy
	244
	245	v_p(k,-1,i) = v(k,-1,i) + dt_3d * c_v * &
	246	( v(k,-1,i) - v_m_s(k,0,i) ) * ddy
	247
	248	w_p(k,-1,i) = w(k,-1,i) + dt_3d * c_w * &
	249	( w(k,-1,i) - w(k,0,i) ) * ddy
	250
	251	!
	252	!-- Save old timelevels for the next timestep
	253	u_m_s(k,:,i) = u(k,-1:1,i)
	254	v_m_s(k,:,i) = v(k,-1:1,i)
	255	w_m_s(k,:,i) = w(k,-1:1,i)
	256
	257	ENDDO
	258	ENDDO
	259
	260	!
	261	!-- Bottom boundary at the outflow
	262	IF ( ibc_uv_b == 0 ) THEN
	263	u_p(nzb,-1,:) = -u_p(nzb+1,-1,:)
	264	v_p(nzb,-1,:) = -v_p(nzb+1,-1,:)
	265	ELSE
	266	u_p(nzb,-1,:) = u_p(nzb+1,-1,:)
	267	v_p(nzb,-1,:) = v_p(nzb+1,-1,:)
[73]	268	ENDIF
[75]	269	w_p(nzb,ny+1,:) = 0.0
[73]	270
[75]	271	!
	272	!-- Top boundary at the outflow
	273	IF ( ibc_uv_t == 0 ) THEN
	274	u_p(nzt+1,-1,:) = ug(nzt+1)
	275	v_p(nzt+1,-1,:) = vg(nzt+1)
	276	ELSE
	277	u_p(nzt+1,-1,:) = u(nzt,-1,:)
	278	v_p(nzt+1,-1,:) = v(nzt,-1,:)
	279	ENDIF
	280	w_p(nzt:nzt+1,-1,:) = 0.0
[73]	281
[75]	282	ENDIF
[73]	283
[75]	284	IF ( outflow_n .AND. &
	285	intermediate_timestep_count == intermediate_timestep_count_max ) &
	286	THEN
[73]	287
[75]	288	c_max = dy / dt_3d
	289
	290	DO i = nxl-1, nxr+1
	291	DO k = nzb+1, nzt+1
	292
[1]	293	!
[75]	294	!-- First calculate the phase speeds for u,v, and w
	295	denom = u_m_n(k,ny,i) - u_m_n(k,ny-1,i)
[73]	296
[75]	297	IF ( denom /= 0.0 ) THEN
	298	c_u = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) / denom
	299	IF ( c_u < 0.0 ) THEN
	300	c_u = 0.0
	301	ELSEIF ( c_u > c_max ) THEN
[73]	302	c_u = c_max
	303	ENDIF
[75]	304	ELSE
	305	c_u = c_max
	306	ENDIF
[73]	307
[75]	308	denom = v_m_n(k,ny,i) - v_m_n(k,ny-1,i)
[73]	309
[75]	310	IF ( denom /= 0.0 ) THEN
	311	c_v = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / denom
	312	IF ( c_v < 0.0 ) THEN
	313	c_v = 0.0
	314	ELSEIF ( c_v > c_max ) THEN
[73]	315	c_v = c_max
	316	ENDIF
[75]	317	ELSE
	318	c_v = c_max
	319	ENDIF
[73]	320
[75]	321	denom = w_m_n(k,ny,i) - w_m_n(k,ny-1,i)
[73]	322
[75]	323	IF ( denom /= 0.0 ) THEN
	324	c_w = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / denom
	325	IF ( c_w < 0.0 ) THEN
	326	c_w = 0.0
	327	ELSEIF ( c_w > c_max ) THEN
[73]	328	c_w = c_max
	329	ENDIF
[75]	330	ELSE
	331	c_w = c_max
	332	ENDIF
[73]	333
	334	!
[75]	335	!-- Calculate the new velocities
	336	u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * c_u * &
	337	( u(k,ny+1,i) - u(k,ny,i) ) * ddy
[73]	338
[75]	339	v_p(k,ny+1,i) = v(k,ny+1,i) - dt_3d * c_v * &
	340	( v(k,ny+1,i) - v(k,ny,i) ) * ddy
[73]	341
[75]	342	w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * c_w * &
	343	( w(k,ny+1,i) - w(k,ny,i) ) * ddy
[73]	344
	345	!
[75]	346	!-- Swap timelevels for the next timestep
	347	u_m_n(k,:,i) = u(k,ny-1:ny+1,i)
	348	v_m_n(k,:,i) = v(k,ny-1:ny+1,i)
	349	w_m_n(k,:,i) = w(k,ny-1:ny+1,i)
[73]	350
[1]	351	ENDDO
[75]	352	ENDDO
[1]	353
	354	!
[75]	355	!-- Bottom boundary at the outflow
	356	IF ( ibc_uv_b == 0 ) THEN
	357	u_p(nzb,ny+1,:) = -u_p(nzb+1,ny+1,:)
	358	v_p(nzb,ny+1,:) = -v_p(nzb+1,ny+1,:)
	359	ELSE
	360	u_p(nzb,ny+1,:) = u_p(nzb+1,ny+1,:)
	361	v_p(nzb,ny+1,:) = v_p(nzb+1,ny+1,:)
	362	ENDIF
	363	w_p(nzb,ny+1,:) = 0.0
[73]	364
	365	!
[75]	366	!-- Top boundary at the outflow
	367	IF ( ibc_uv_t == 0 ) THEN
	368	u_p(nzt+1,ny+1,:) = ug(nzt+1)
	369	v_p(nzt+1,ny+1,:) = vg(nzt+1)
	370	ELSE
	371	u_p(nzt+1,ny+1,:) = u_p(nzt,nyn+1,:)
	372	v_p(nzt+1,ny+1,:) = v_p(nzt,nyn+1,:)
[1]	373	ENDIF
[75]	374	w_p(nzt:nzt+1,ny+1,:) = 0.0
[1]	375
[75]	376	ENDIF
	377
	378	IF ( outflow_l .AND. &
	379	intermediate_timestep_count == intermediate_timestep_count_max ) &
	380	THEN
	381
	382	c_max = dx / dt_3d
	383
	384	DO j = nys-1, nyn+1
	385	DO k = nzb+1, nzt+1
	386
[1]	387	!
[75]	388	!-- First calculate the phase speeds for u,v, and w
	389	denom = u_m_l(k,j,0) - u_m_l(k,j,1)
	390
	391	IF ( denom /= 0.0 ) THEN
	392	c_u = -c_max * ( u(k,j,0) - u_m_r(k,j,0) ) / denom
	393	IF ( c_u > 0.0 ) THEN
	394	c_u = 0.0
	395	ELSEIF ( c_u < -c_max ) THEN
	396	c_u = -c_max
	397	ENDIF
	398	ELSE
	399	c_u = -c_max
	400	ENDIF
	401
	402	denom = v_m_l(k,j,0) - v_m_l(k,j,1)
	403
	404	IF ( denom /= 0.0 ) THEN
	405	c_v = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) / denom
	406	IF ( c_v < 0.0 ) THEN
	407	c_v = 0.0
	408	ELSEIF ( c_v > c_max ) THEN
	409	c_v = c_max
	410	ENDIF
	411	ELSE
	412	c_v = c_max
	413	ENDIF
	414
	415	denom = w_m_l(k,j,0) - w_m_l(k,j,1)
	416
	417	IF ( denom /= 0.0 ) THEN
	418	c_w = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) / denom
	419	IF ( c_w < 0.0 ) THEN
	420	c_w = 0.0
	421	ELSEIF ( c_w > c_max ) THEN
	422	c_w = c_max
	423	ENDIF
	424	ELSE
	425	c_w = c_max
	426	ENDIF
	427
[73]	428	!
[75]	429	!-- Calculate the new velocities
	430	u_p(k,j,-1) = u(k,j,-1) + dt_3d * c_u * &
	431	( u(k,j,-1) - u(k,j,0) ) * ddx
	432
	433	v_p(k,j,-1) = v(k,j,-1) + dt_3d * c_v * &
	434	( v(k,j,-1) - v(k,j,0) ) * ddx
	435
	436	w_p(k,j,-1) = w(k,j,-1) + dt_3d * c_w * &
	437	( w(k,j,-1) - w(k,j,0) ) * ddx
	438
	439	!
	440	!-- Swap timelevels for the next timestep
	441	u_m_l(k,j,:) = u(k,j,-1:1)
	442	v_m_l(k,j,:) = v(k,j,-1:1)
	443	w_m_l(k,j,:) = w(k,j,-1:1)
	444
	445	ENDDO
	446	ENDDO
	447
	448	!
	449	!-- Bottom boundary at the outflow
	450	IF ( ibc_uv_b == 0 ) THEN
	451	u_p(nzb,:,-1) = -u_p(nzb+1,:,-1)
	452	v_p(nzb,:,-1) = -v_p(nzb+1,:,-1)
	453	ELSE
	454	u_p(nzb,:,-1) = u_p(nzb+1,:,-1)
	455	v_p(nzb,:,-1) = v_p(nzb+1,:,-1)
[1]	456	ENDIF
[75]	457	w_p(nzb,:,-1) = 0.0
[1]	458
[75]	459	!
	460	!-- Top boundary at the outflow
	461	IF ( ibc_uv_t == 0 ) THEN
	462	u_p(nzt+1,:,-1) = ug(nzt+1)
	463	v_p(nzt+1,:,-1) = vg(nzt+1)
	464	ELSE
	465	u_p(nzt+1,:,-1) = u_p(nzt,:,-1)
	466	v_p(nzt+1,:,-1) = v_p(nzt,:,-1)
	467	ENDIF
	468	w_p(nzt:nzt+1,:,-1) = 0.0
[73]	469
[75]	470	ENDIF
[73]	471
[75]	472	IF ( outflow_r .AND. &
	473	intermediate_timestep_count == intermediate_timestep_count_max ) &
	474	THEN
[73]	475
[75]	476	c_max = dx / dt_3d
	477
	478	DO j = nys-1, nyn+1
	479	DO k = nzb+1, nzt+1
	480
[1]	481	!
[75]	482	!-- First calculate the phase speeds for u,v, and w
	483	denom = u_m_r(k,j,nx) - u_m_r(k,j,nx-1)
[73]	484
[75]	485	IF ( denom /= 0.0 ) THEN
	486	c_u = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) ) / denom
	487	IF ( c_u < 0.0 ) THEN
	488	c_u = 0.0
	489	ELSEIF ( c_u > c_max ) THEN
[73]	490	c_u = c_max
	491	ENDIF
[75]	492	ELSE
	493	c_u = c_max
	494	ENDIF
[73]	495
[75]	496	denom = v_m_r(k,j,nx) - v_m_r(k,j,nx-1)
[73]	497
[75]	498	IF ( denom /= 0.0 ) THEN
	499	c_v = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) / denom
	500	IF ( c_v < 0.0 ) THEN
	501	c_v = 0.0
	502	ELSEIF ( c_v > c_max ) THEN
[73]	503	c_v = c_max
	504	ENDIF
[75]	505	ELSE
	506	c_v = c_max
	507	ENDIF
[73]	508
[75]	509	denom = w_m_r(k,j,nx) - w_m_r(k,j,nx-1)
[73]	510
[75]	511	IF ( denom /= 0.0 ) THEN
	512	c_w = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) / denom
	513	IF ( c_w < 0.0 ) THEN
	514	c_w = 0.0
	515	ELSEIF ( c_w > c_max ) THEN
[73]	516	c_w = c_max
	517	ENDIF
[75]	518	ELSE
	519	c_w = c_max
	520	ENDIF
[73]	521
	522	!
[75]	523	!-- Calculate the new velocities
	524	u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * c_u * &
	525	( u(k,j,nx+1) - u(k,j,nx) ) * ddx
[73]	526
[75]	527	v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * c_v * &
	528	( v(k,j,nx+1) - v(k,j,nx) ) * ddx
[73]	529
[75]	530	w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * c_w * &
	531	( w(k,j,nx+1) - w(k,j,nx) ) * ddx
[73]	532
	533	!
[75]	534	!-- Swap timelevels for the next timestep
	535	u_m_r(k,j,:) = u(k,j,nx-1:nx+1)
	536	v_m_r(k,j,:) = v(k,j,nx-1:nx+1)
	537	w_m_r(k,j,:) = w(k,j,nx-1:nx+1)
[73]	538
	539	ENDDO
[75]	540	ENDDO
[73]	541
	542	!
[75]	543	!-- Bottom boundary at the outflow
	544	IF ( ibc_uv_b == 0 ) THEN
	545	u_p(nzb,:,nx+1) = -u_p(nzb+1,:,nx+1)
	546	v_p(nzb,:,nx+1) = -v_p(nzb+1,:,nx+1)
	547	ELSE
	548	u_p(nzb,:,nx+1) = u_p(nzb+1,:,nx+1)
	549	v_p(nzb,:,nx+1) = v_p(nzb+1,:,nx+1)
	550	ENDIF
	551	w_p(nzb,:,nx+1) = 0.0
[73]	552
	553	!
[75]	554	!-- Top boundary at the outflow
	555	IF ( ibc_uv_t == 0 ) THEN
	556	u_p(nzt+1,:,nx+1) = ug(nzt+1)
	557	v_p(nzt+1,:,nx+1) = vg(nzt+1)
	558	ELSE
	559	u_p(nzt+1,:,nx+1) = u_p(nzt,:,nx+1)
	560	v_p(nzt+1,:,nx+1) = v_p(nzt,:,nx+1)
[1]	561	ENDIF
[75]	562	w(nzt:nzt+1,:,nx+1) = 0.0
[1]	563
	564	ENDIF
	565
	566
	567	END SUBROUTINE boundary_conds

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

Context Navigation

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

Download in other formats: