Home

Context Navigation

init_3d_model.f90 @ 141

Last change on this file since 141 was 139, checked in by raasch, 17 years ago

New:
---

Plant canopy model of Watanabe (2004,BLM 112,307-341) added.
It can be switched on by the inipar parameter plant_canopy.
The inipar parameter canopy_mode can be used to prescribe a
plant canopy type. The default case is a homogeneous plant
canopy. Heterogeneous distributions of the leaf area
density and the canopy drag coefficient can be defined in the
new routine user_init_plant_canopy (user_interface).
The inipar parameters lad_surface, lad_vertical_gradient and
lad_vertical_gradient_level can be used in order to
prescribe the vertical profile of leaf area density. The
inipar parameter drag_coefficient determines the canopy
drag coefficient.
Finally, the inipar parameter pch_index determines the
index of the upper boundary of the plant canopy.

Allow new case bc_uv_t = 'dirichlet_0' for channel flow.

For unknown variables (CASE DEFAULT) call new subroutine user_data_output_dvrp

Pressure boundary conditions for vertical walls added to the multigrid solver.
They are applied using new wall flag arrays (wall_flags_..) which are defined
for each grid level. New argument gls added to routine user_init_grid
(user_interface).

Frequence of sorting particles can be controlled with new particles_par
parameter dt_sort_particles. Sorting is moved from the SGS timestep loop in
advec_particles after the end of this loop.

advec_particles, check_parameters, data_output_dvrp, header, init_3d_model, init_grid, init_particles, init_pegrid, modules, package_parin, parin, plant_canopy_model, read_var_list, read_3d_binary, user_interface, write_var_list, write_3d_binary

Changed:

Redefine initial nzb_local as the actual total size of topography (later the
extent of topography in nzb_local is reduced by 1dx at the E topography walls
and by 1dy at the N topography walls to form the basis for nzb_s_inner);
for consistency redefine 'single_building' case.

Vertical profiles now based on nzb_s_inner; they are divided by
ngp_2dh_s_inner (scalars, procucts of scalars) and ngp_2dh (staggered velocity
components and their products, procucts of scalars and velocity components),
respectively.

Allow two instead of one digit to specify isosurface and slicer variables.

Status of 3D-volume NetCDF data file only depends on switch netcdf_64bit_3d (check_open)

prognostic_equations include the respective wall_*flux in the parameter list of
calls of diffusion_s. Same as before, only the values of wall_heatflux(0:4)
can be assigned. At present, wall_humidityflux, wall_qflux, wall_salinityflux,
and wall_scalarflux are kept zero. diffusion_s uses the respective wall_*flux
instead of wall_heatflux. This update serves two purposes:

it avoids errors in calculations with humidity/scalar/salinity and prescribed

non-zero wall_heatflux,

it prepares PALM for a possible assignment of wall fluxes of

humidity/scalar/salinity in a future release.

buoyancy, check_open, data_output_dvrp, diffusion_s, diffusivities, flow_statistics, header, init_3d_model, init_dvrp, init_grid, modules, prognostic_equations

Errors:

Bugfix: summation of sums_l_l in diffusivities.

Several bugfixes in the ocean part: Initial density rho is calculated
(init_ocean). Error in initializing u_init and v_init removed
(check_parameters). Calculation of density flux now starts from
nzb+1 (production_e).

Bugfix: pleft/pright changed to pnorth/psouth in sendrecv of particle tail
numbers along y, small bugfixes in the SGS part (advec_particles)

Bugfix: model_string needed a default value (combine_plot_fields)

Bugfix: wavenumber calculation for even nx in routines maketri (poisfft)

Bugfix: assignment of fluxes at walls

Bugfix: absolute value of f must be used when calculating the Blackadar mixing length (init_1d_model)

advec_particles, check_parameters, combine_plot_fields, diffusion_s, diffusivities, init_ocean, init_1d_model, poisfft, production_e

Property svn:keywords set to Id

File size: 40.8 KB

Rev	Line
[1]	1	#if defined( __ibmy_special )
	2	@PROCESS NOOPTimize
	3	#endif
	4	SUBROUTINE init_3d_model
	5
	6	!------------------------------------------------------------------------------!
	7	! Actual revisions:
	8	! -----------------
[139]	9	!
	10	!
	11	! Former revisions:
	12	! -----------------
	13	! $Id: init_3d_model.f90 139 2007-11-29 09:37:41Z raasch $
	14	!
	15	! 138 2007-11-28 10:03:58Z letzel
[132]	16	! New counter ngp_2dh_s_inner.
	17	! Allow new case bc_uv_t = 'dirichlet_0' for channel flow.
	18	! Corrected calculation of initial volume flow for 'set_1d-model_profiles' and
	19	! 'set_constant_profiles' in case of buildings in the reference cross-sections.
[77]	20	!
[110]	21	! 108 2007-08-24 15:10:38Z letzel
	22	! Flux initialization in case of coupled runs, +momentum fluxes at top boundary,
	23	! +arrays for phase speed c_u, c_v, c_w, indices for u\|v\|w_m_l\|r changed
	24	! +qswst_remote in case of atmosphere model with humidity coupled to ocean
	25	! Rayleigh damping for ocean, optionally calculate km and kh from initial
	26	! TKE e_init
	27	!
[98]	28	! 97 2007-06-21 08:23:15Z raasch
	29	! Initialization of salinity, call of init_ocean
	30	!
[90]	31	! 87 2007-05-22 15:46:47Z raasch
	32	! var_hom and var_sum renamed pr_palm
	33	!
[77]	34	! 75 2007-03-22 09:54:05Z raasch
[73]	35	! Arrays for radiation boundary conditions are allocated (u_m_l, u_m_r, etc.),
	36	! bugfix for cases with the outflow damping layer extending over more than one
[75]	37	! subdomain, moisture renamed humidity,
	38	! new initializing action "by_user" calls user_init_3d_model,
[72]	39	! precipitation_amount/rate, ts_value are allocated, +module netcdf_control,
[51]	40	! initial velocities at nzb+1 are regarded for volume
	41	! flow control in case they have been set zero before (to avoid small timesteps)
[75]	42	! -uvmean_outflow, uxrp, vynp eliminated
[1]	43	!
[39]	44	! 19 2007-02-23 04:53:48Z raasch
	45	! +handling of top fluxes
	46	!
[3]	47	! RCS Log replace by Id keyword, revision history cleaned up
	48	!
[1]	49	! Revision 1.49 2006/08/22 15:59:07 raasch
	50	! No optimization of this file on the ibmy (Yonsei Univ.)
	51	!
	52	! Revision 1.1 1998/03/09 16:22:22 raasch
	53	! Initial revision
	54	!
	55	!
	56	! Description:
	57	! ------------
	58	! Allocation of arrays and initialization of the 3D model via
	59	! a) pre-run the 1D model
	60	! or
	61	! b) pre-set constant linear profiles
	62	! or
	63	! c) read values of a previous run
	64	!------------------------------------------------------------------------------!
	65
	66	USE arrays_3d
	67	USE averaging
[72]	68	USE cloud_parameters
[1]	69	USE constants
	70	USE control_parameters
	71	USE cpulog
	72	USE indices
	73	USE interfaces
	74	USE model_1d
[51]	75	USE netcdf_control
[1]	76	USE particle_attributes
	77	USE pegrid
	78	USE profil_parameter
	79	USE random_function_mod
	80	USE statistics
	81
	82	IMPLICIT NONE
	83
	84	INTEGER :: i, j, k, sr
	85
	86	INTEGER, DIMENSION(:), ALLOCATABLE :: ngp_2dh_l, ngp_3d_inner_l
	87
[132]	88	INTEGER, DIMENSION(:,:), ALLOCATABLE :: ngp_2dh_outer_l, &
	89	ngp_2dh_s_inner_l
[1]	90
	91	REAL, DIMENSION(1:2) :: volume_flow_area_l, volume_flow_initial_l
	92
	93
	94	!
	95	!-- Allocate arrays
	96	ALLOCATE( ngp_2dh(0:statistic_regions), ngp_2dh_l(0:statistic_regions), &
	97	ngp_3d(0:statistic_regions), &
	98	ngp_3d_inner(0:statistic_regions), &
	99	ngp_3d_inner_l(0:statistic_regions), &
	100	sums_divnew_l(0:statistic_regions), &
	101	sums_divold_l(0:statistic_regions) )
[75]	102	ALLOCATE( rdf(nzb+1:nzt) )
[87]	103	ALLOCATE( hom_sum(nzb:nzt+1,pr_palm+max_pr_user,0:statistic_regions), &
[1]	104	ngp_2dh_outer(nzb:nzt+1,0:statistic_regions), &
	105	ngp_2dh_outer_l(nzb:nzt+1,0:statistic_regions), &
[132]	106	ngp_2dh_s_inner(nzb:nzt+1,0:statistic_regions), &
	107	ngp_2dh_s_inner_l(nzb:nzt+1,0:statistic_regions), &
[1]	108	rmask(nys-1:nyn+1,nxl-1:nxr+1,0:statistic_regions), &
[87]	109	sums(nzb:nzt+1,pr_palm+max_pr_user), &
	110	sums_l(nzb:nzt+1,pr_palm+max_pr_user,0:threads_per_task-1), &
[1]	111	sums_l_l(nzb:nzt+1,0:statistic_regions,0:threads_per_task-1), &
	112	sums_up_fraction_l(10,3,0:statistic_regions), &
[48]	113	sums_wsts_bc_l(nzb:nzt+1,0:statistic_regions), &
	114	ts_value(var_ts,0:statistic_regions) )
[1]	115	ALLOCATE( km_damp_x(nxl-1:nxr+1), km_damp_y(nys-1:nyn+1) )
	116
[19]	117	ALLOCATE( rif_1(nys-1:nyn+1,nxl-1:nxr+1), shf_1(nys-1:nyn+1,nxl-1:nxr+1), &
	118	ts(nys-1:nyn+1,nxl-1:nxr+1), tswst_1(nys-1:nyn+1,nxl-1:nxr+1), &
	119	us(nys-1:nyn+1,nxl-1:nxr+1), usws_1(nys-1:nyn+1,nxl-1:nxr+1), &
[102]	120	uswst_1(nys-1:nyn+1,nxl-1:nxr+1), &
	121	vsws_1(nys-1:nyn+1,nxl-1:nxr+1), &
	122	vswst_1(nys-1:nyn+1,nxl-1:nxr+1), z0(nys-1:nyn+1,nxl-1:nxr+1) )
[1]	123
	124	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	125	!
	126	!-- Leapfrog scheme needs two timelevels of diffusion quantities
[19]	127	ALLOCATE( rif_2(nys-1:nyn+1,nxl-1:nxr+1), &
	128	shf_2(nys-1:nyn+1,nxl-1:nxr+1), &
	129	tswst_2(nys-1:nyn+1,nxl-1:nxr+1), &
	130	usws_2(nys-1:nyn+1,nxl-1:nxr+1), &
[102]	131	uswst_2(nys-1:nyn+1,nxl-1:nxr+1), &
	132	vswst_2(nys-1:nyn+1,nxl-1:nxr+1), &
[1]	133	vsws_2(nys-1:nyn+1,nxl-1:nxr+1) )
	134	ENDIF
	135
[75]	136	ALLOCATE( d(nzb+1:nzta,nys:nyna,nxl:nxra), &
	137	e_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
	138	e_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
	139	e_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
	140	kh_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
	141	km_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
	142	p(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
	143	pt_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
	144	pt_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
	145	pt_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
	146	tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
	147	u_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
	148	u_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
	149	u_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
	150	v_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
	151	v_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
	152	v_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
	153	w_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
	154	w_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
[1]	155	w_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
	156
	157	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	158	ALLOCATE( kh_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
	159	km_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
	160	ENDIF
	161
[75]	162	IF ( humidity .OR. passive_scalar ) THEN
[1]	163	!
[75]	164	!-- 2D-humidity/scalar arrays
[1]	165	ALLOCATE ( qs(nys-1:nyn+1,nxl-1:nxr+1), &
[19]	166	qsws_1(nys-1:nyn+1,nxl-1:nxr+1), &
	167	qswst_1(nys-1:nyn+1,nxl-1:nxr+1) )
[1]	168
	169	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
[19]	170	ALLOCATE( qsws_2(nys-1:nyn+1,nxl-1:nxr+1), &
	171	qswst_2(nys-1:nyn+1,nxl-1:nxr+1) )
[1]	172	ENDIF
	173	!
[75]	174	!-- 3D-humidity/scalar arrays
[1]	175	ALLOCATE( q_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
	176	q_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
	177	q_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
	178
	179	!
[75]	180	!-- 3D-arrays needed for humidity only
	181	IF ( humidity ) THEN
[1]	182	ALLOCATE( vpt_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
	183
	184	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	185	ALLOCATE( vpt_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
	186	ENDIF
	187
	188	IF ( cloud_physics ) THEN
	189	!
	190	!-- Liquid water content
	191	ALLOCATE ( ql_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
[72]	192	!
	193	!-- Precipitation amount and rate (only needed if output is switched)
	194	ALLOCATE( precipitation_amount(nys-1:nyn+1,nxl-1:nxr+1), &
	195	precipitation_rate(nys-1:nyn+1,nxl-1:nxr+1) )
[1]	196	ENDIF
	197
	198	IF ( cloud_droplets ) THEN
	199	!
	200	!-- Liquid water content, change in liquid water content,
	201	!-- real volume of particles (with weighting), volume of particles
	202	ALLOCATE ( ql_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
	203	ql_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
	204	ql_v(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
	205	ql_vp(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
	206	ENDIF
	207
	208	ENDIF
	209
	210	ENDIF
	211
[94]	212	IF ( ocean ) THEN
[95]	213	ALLOCATE( saswsb_1(nys-1:nyn+1,nxl-1:nxr+1), &
	214	saswst_1(nys-1:nyn+1,nxl-1:nxr+1) )
[96]	215	ALLOCATE( rho_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
	216	sa_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
	217	sa_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
[94]	218	sa_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
[96]	219	rho => rho_1 ! routine calc_mean_profile requires density to be a
	220	! pointer
[108]	221	IF ( humidity_remote ) THEN
	222	ALLOCATE( qswst_remote(nys-1:nyn+1,nxl-1:nxr+1) )
	223	qswst_remote = 0.0
	224	ENDIF
[94]	225	ENDIF
	226
[1]	227	!
	228	!-- 3D-array for storing the dissipation, needed for calculating the sgs
	229	!-- particle velocities
	230	IF ( use_sgs_for_particles ) THEN
	231	ALLOCATE ( diss(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
	232	ENDIF
	233
	234	IF ( dt_dosp /= 9999999.9 ) THEN
	235	ALLOCATE( spectrum_x( 1:nx/2, 1:10, 1:10 ), &
	236	spectrum_y( 1:ny/2, 1:10, 1:10 ) )
	237	ENDIF
	238
	239	!
[138]	240	!-- 3D-arrays for the leaf area density and the canopy drag coefficient
	241	IF ( plant_canopy ) THEN
	242	ALLOCATE ( lad_s(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
	243	lad_u(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
	244	lad_v(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
	245	lad_w(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
	246	cdc(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
	247	ENDIF
	248
	249	!
[51]	250	!-- 4D-array for storing the Rif-values at vertical walls
	251	IF ( topography /= 'flat' ) THEN
	252	ALLOCATE( rif_wall(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1,1:4) )
	253	rif_wall = 0.0
	254	ENDIF
	255
	256	!
	257	!-- Velocities at nzb+1 needed for volume flow control
	258	IF ( conserve_volume_flow ) THEN
	259	ALLOCATE( u_nzb_p1_for_vfc(nys:nyn), v_nzb_p1_for_vfc(nxl:nxr) )
	260	u_nzb_p1_for_vfc = 0.0
	261	v_nzb_p1_for_vfc = 0.0
	262	ENDIF
	263
	264	!
[106]	265	!-- Arrays to store velocity data from t-dt and the phase speeds which
	266	!-- are needed for radiation boundary conditions
[73]	267	IF ( outflow_l ) THEN
[106]	268	ALLOCATE( u_m_l(nzb:nzt+1,nys-1:nyn+1,1:2), &
	269	v_m_l(nzb:nzt+1,nys-1:nyn+1,0:1), &
	270	w_m_l(nzb:nzt+1,nys-1:nyn+1,0:1) )
[73]	271	ENDIF
	272	IF ( outflow_r ) THEN
[106]	273	ALLOCATE( u_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx), &
	274	v_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx), &
	275	w_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx) )
[73]	276	ENDIF
[106]	277	IF ( outflow_l .OR. outflow_r ) THEN
	278	ALLOCATE( c_u(nzb:nzt+1,nys-1:nyn+1), c_v(nzb:nzt+1,nys-1:nyn+1), &
	279	c_w(nzb:nzt+1,nys-1:nyn+1) )
	280	ENDIF
[73]	281	IF ( outflow_s ) THEN
[106]	282	ALLOCATE( u_m_s(nzb:nzt+1,0:1,nxl-1:nxr+1), &
	283	v_m_s(nzb:nzt+1,1:2,nxl-1:nxr+1), &
	284	w_m_s(nzb:nzt+1,0:1,nxl-1:nxr+1) )
[73]	285	ENDIF
	286	IF ( outflow_n ) THEN
[106]	287	ALLOCATE( u_m_n(nzb:nzt+1,ny-1:ny,nxl-1:nxr+1), &
	288	v_m_n(nzb:nzt+1,ny-1:ny,nxl-1:nxr+1), &
	289	w_m_n(nzb:nzt+1,ny-1:ny,nxl-1:nxr+1) )
[73]	290	ENDIF
[106]	291	IF ( outflow_s .OR. outflow_n ) THEN
	292	ALLOCATE( c_u(nzb:nzt+1,nxl-1:nxr+1), c_v(nzb:nzt+1,nxl-1:nxr+1), &
	293	c_w(nzb:nzt+1,nxl-1:nxr+1) )
	294	ENDIF
[73]	295
	296	!
[1]	297	!-- Initial assignment of the pointers
	298	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	299
[19]	300	rif_m => rif_1; rif => rif_2
	301	shf_m => shf_1; shf => shf_2
	302	tswst_m => tswst_1; tswst => tswst_2
	303	usws_m => usws_1; usws => usws_2
[102]	304	uswst_m => uswst_1; uswst => uswst_2
[19]	305	vsws_m => vsws_1; vsws => vsws_2
[102]	306	vswst_m => vswst_1; vswst => vswst_2
[1]	307	e_m => e_1; e => e_2; e_p => e_3; te_m => e_3
	308	kh_m => kh_1; kh => kh_2
	309	km_m => km_1; km => km_2
	310	pt_m => pt_1; pt => pt_2; pt_p => pt_3; tpt_m => pt_3
	311	u_m => u_1; u => u_2; u_p => u_3; tu_m => u_3
	312	v_m => v_1; v => v_2; v_p => v_3; tv_m => v_3
	313	w_m => w_1; w => w_2; w_p => w_3; tw_m => w_3
	314
[75]	315	IF ( humidity .OR. passive_scalar ) THEN
[19]	316	qsws_m => qsws_1; qsws => qsws_2
	317	qswst_m => qswst_1; qswst => qswst_2
[1]	318	q_m => q_1; q => q_2; q_p => q_3; tq_m => q_3
[75]	319	IF ( humidity ) vpt_m => vpt_1; vpt => vpt_2
[1]	320	IF ( cloud_physics ) ql => ql_1
	321	IF ( cloud_droplets ) THEN
	322	ql => ql_1
	323	ql_c => ql_2
	324	ENDIF
	325	ENDIF
	326
	327	ELSE
	328
[19]	329	rif => rif_1
	330	shf => shf_1
	331	tswst => tswst_1
	332	usws => usws_1
[102]	333	uswst => uswst_1
[19]	334	vsws => vsws_1
[102]	335	vswst => vswst_1
[19]	336	e => e_1; e_p => e_2; te_m => e_3; e_m => e_3
	337	kh => kh_1
	338	km => km_1
	339	pt => pt_1; pt_p => pt_2; tpt_m => pt_3; pt_m => pt_3
	340	u => u_1; u_p => u_2; tu_m => u_3; u_m => u_3
	341	v => v_1; v_p => v_2; tv_m => v_3; v_m => v_3
	342	w => w_1; w_p => w_2; tw_m => w_3; w_m => w_3
[1]	343
[75]	344	IF ( humidity .OR. passive_scalar ) THEN
[1]	345	qsws => qsws_1
[19]	346	qswst => qswst_1
[94]	347	q => q_1; q_p => q_2; tq_m => q_3; q_m => q_3
[75]	348	IF ( humidity ) vpt => vpt_1
[1]	349	IF ( cloud_physics ) ql => ql_1
	350	IF ( cloud_droplets ) THEN
	351	ql => ql_1
	352	ql_c => ql_2
	353	ENDIF
	354	ENDIF
	355
[94]	356	IF ( ocean ) THEN
[95]	357	saswsb => saswsb_1
[94]	358	saswst => saswst_1
	359	sa => sa_1; sa_p => sa_2; tsa_m => sa_3
	360	ENDIF
	361
[1]	362	ENDIF
	363
	364	!
	365	!-- Initialize model variables
	366	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
	367	!
	368	!-- First model run of a possible job queue.
	369	!-- Initial profiles of the variables must be computes.
	370	IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 ) THEN
	371	!
	372	!-- Use solutions of the 1D model as initial profiles,
	373	!-- start 1D model
	374	CALL init_1d_model
	375	!
	376	!-- Transfer initial profiles to the arrays of the 3D model
	377	DO i = nxl-1, nxr+1
	378	DO j = nys-1, nyn+1
	379	e(:,j,i) = e1d
	380	kh(:,j,i) = kh1d
	381	km(:,j,i) = km1d
	382	pt(:,j,i) = pt_init
	383	u(:,j,i) = u1d
	384	v(:,j,i) = v1d
	385	ENDDO
	386	ENDDO
	387
[75]	388	IF ( humidity .OR. passive_scalar ) THEN
[1]	389	DO i = nxl-1, nxr+1
	390	DO j = nys-1, nyn+1
	391	q(:,j,i) = q_init
	392	ENDDO
	393	ENDDO
	394	ENDIF
	395
	396	IF ( .NOT. constant_diffusion ) THEN
	397	DO i = nxl-1, nxr+1
	398	DO j = nys-1, nyn+1
	399	e(:,j,i) = e1d
	400	ENDDO
	401	ENDDO
	402	!
	403	!-- Store initial profiles for output purposes etc.
	404	hom(:,1,25,:) = SPREAD( l1d, 2, statistic_regions+1 )
	405
	406	IF ( prandtl_layer ) THEN
	407	rif = rif1d(nzb+1)
	408	ts = 0.0 ! could actually be computed more accurately in the
	409	! 1D model. Update when opportunity arises.
	410	us = us1d
	411	usws = usws1d
	412	vsws = vsws1d
	413	ELSE
	414	ts = 0.0 ! must be set, because used in
	415	rif = 0.0 ! flowste
	416	us = 0.0
	417	usws = 0.0
	418	vsws = 0.0
	419	ENDIF
	420
	421	ELSE
	422	e = 0.0 ! must be set, because used in
	423	rif = 0.0 ! flowste
	424	ts = 0.0
	425	us = 0.0
	426	usws = 0.0
	427	vsws = 0.0
	428	ENDIF
[102]	429	uswst = top_momentumflux_u
	430	vswst = top_momentumflux_v
[1]	431
	432	!
	433	!-- In every case qs = 0.0 (see also pt)
	434	!-- This could actually be computed more accurately in the 1D model.
	435	!-- Update when opportunity arises!
[75]	436	IF ( humidity .OR. passive_scalar ) qs = 0.0
[1]	437
	438	!
	439	!-- inside buildings set velocities back to zero
	440	IF ( topography /= 'flat' ) THEN
	441	DO i = nxl-1, nxr+1
	442	DO j = nys-1, nyn+1
	443	u(nzb:nzb_u_inner(j,i),j,i) = 0.0
	444	v(nzb:nzb_v_inner(j,i),j,i) = 0.0
	445	ENDDO
	446	ENDDO
[132]	447	IF ( conserve_volume_flow ) THEN
	448	IF ( nxr == nx ) THEN
	449	DO j = nys, nyn
	450	DO k = nzb + 1, nzb_u_inner(j,nx)
	451	u_nzb_p1_for_vfc(j) = u1d(k) * dzu(k)
	452	ENDDO
	453	ENDDO
	454	ENDIF
	455	IF ( nyn == ny ) THEN
	456	DO i = nxl, nxr
	457	DO k = nzb + 1, nzb_v_inner(ny,i)
	458	v_nzb_p1_for_vfc(i) = v1d(k) * dzu(k)
	459	ENDDO
	460	ENDDO
	461	ENDIF
	462	ENDIF
[1]	463	!
	464	!-- WARNING: The extra boundary conditions set after running the
	465	!-- ------- 1D model impose an error on the divergence one layer
	466	!-- below the topography; need to correct later
	467	!-- ATTENTION: Provisional correction for Piacsek & Williams
	468	!-- --------- advection scheme: keep u and v zero one layer below
	469	!-- the topography.
	470	IF ( ibc_uv_b == 0 ) THEN
	471	!
	472	!-- Satisfying the Dirichlet condition with an extra layer below
	473	!-- the surface where the u and v component change their sign.
	474	DO i = nxl-1, nxr+1
	475	DO j = nys-1, nyn+1
	476	IF ( nzb_u_inner(j,i) == 0 ) u(0,j,i) = -u(1,j,i)
	477	IF ( nzb_v_inner(j,i) == 0 ) v(0,j,i) = -v(1,j,i)
	478	ENDDO
	479	ENDDO
	480
	481	ELSE
	482	!
	483	!-- Neumann condition
	484	DO i = nxl-1, nxr+1
	485	DO j = nys-1, nyn+1
	486	IF ( nzb_u_inner(j,i) == 0 ) u(0,j,i) = u(1,j,i)
	487	IF ( nzb_v_inner(j,i) == 0 ) v(0,j,i) = v(1,j,i)
	488	ENDDO
	489	ENDDO
	490
	491	ENDIF
	492
	493	ENDIF
	494
	495	ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 ) &
	496	THEN
	497	!
	498	!-- Use constructed initial profiles (velocity constant with height,
	499	!-- temperature profile with constant gradient)
	500	DO i = nxl-1, nxr+1
	501	DO j = nys-1, nyn+1
	502	pt(:,j,i) = pt_init
	503	u(:,j,i) = u_init
	504	v(:,j,i) = v_init
	505	ENDDO
	506	ENDDO
[75]	507
[1]	508	!
[51]	509	!-- Set initial horizontal velocities at the lowest computational grid levels
	510	!-- to zero in order to avoid too small time steps caused by the diffusion
[1]	511	!-- limit in the initial phase of a run (at k=1, dz/2 occurs in the
[51]	512	!-- limiting formula!). The original values are stored to be later used for
	513	!-- volume flow control.
[1]	514	DO i = nxl-1, nxr+1
	515	DO j = nys-1, nyn+1
	516	u(nzb:nzb_u_inner(j,i)+1,j,i) = 0.0
	517	v(nzb:nzb_v_inner(j,i)+1,j,i) = 0.0
	518	ENDDO
	519	ENDDO
[51]	520	IF ( conserve_volume_flow ) THEN
	521	IF ( nxr == nx ) THEN
	522	DO j = nys, nyn
[132]	523	DO k = nzb + 1, nzb_u_inner(j,nx) + 1
	524	u_nzb_p1_for_vfc(j) = u_init(k) * dzu(k)
	525	ENDDO
[51]	526	ENDDO
	527	ENDIF
	528	IF ( nyn == ny ) THEN
	529	DO i = nxl, nxr
[132]	530	DO k = nzb + 1, nzb_v_inner(ny,i) + 1
	531	v_nzb_p1_for_vfc(i) = v_init(k) * dzu(k)
	532	ENDDO
[51]	533	ENDDO
	534	ENDIF
	535	ENDIF
[1]	536
[75]	537	IF ( humidity .OR. passive_scalar ) THEN
[1]	538	DO i = nxl-1, nxr+1
	539	DO j = nys-1, nyn+1
	540	q(:,j,i) = q_init
	541	ENDDO
	542	ENDDO
	543	ENDIF
	544
[94]	545	IF ( ocean ) THEN
	546	DO i = nxl-1, nxr+1
	547	DO j = nys-1, nyn+1
	548	sa(:,j,i) = sa_init
	549	ENDDO
	550	ENDDO
	551	ENDIF
[1]	552
	553	IF ( constant_diffusion ) THEN
	554	km = km_constant
	555	kh = km / prandtl_number
[108]	556	e = 0.0
	557	ELSEIF ( e_init > 0.0 ) THEN
	558	DO k = nzb+1, nzt
	559	km(k,:,:) = 0.1 * l_grid(k) * SQRT( e_init )
	560	ENDDO
	561	km(nzb,:,:) = km(nzb+1,:,:)
	562	km(nzt+1,:,:) = km(nzt,:,:)
	563	kh = km / prandtl_number
	564	e = e_init
[1]	565	ELSE
[108]	566	IF ( .NOT. ocean ) THEN
	567	kh = 0.01 ! there must exist an initial diffusion, because
	568	km = 0.01 ! otherwise no TKE would be produced by the
	569	! production terms, as long as not yet
	570	! e = (u/cm)*2 at k=nzb+1
	571	ELSE
	572	kh = 0.00001
	573	km = 0.00001
	574	ENDIF
	575	e = 0.0
[1]	576	ENDIF
[102]	577	rif = 0.0
	578	ts = 0.0
	579	us = 0.0
	580	usws = 0.0
	581	uswst = top_momentumflux_u
	582	vsws = 0.0
	583	vswst = top_momentumflux_v
[75]	584	IF ( humidity .OR. passive_scalar ) qs = 0.0
[1]	585
	586	!
	587	!-- Compute initial temperature field and other constants used in case
	588	!-- of a sloping surface
	589	IF ( sloping_surface ) CALL init_slope
	590
[46]	591	ELSEIF ( INDEX(initializing_actions, 'by_user') /= 0 ) &
	592	THEN
	593	!
	594	!-- Initialization will completely be done by the user
	595	CALL user_init_3d_model
	596
[1]	597	ENDIF
	598
	599	!
[132]	600	!-- apply channel flow boundary condition
	601	IF ( TRIM( bc_uv_t ) == 'dirichlet_0' ) THEN
	602
	603	u(nzt+1,:,:) = 0.0
	604	v(nzt+1,:,:) = 0.0
	605
	606	!-- for the Dirichlet condition to be correctly applied at the top, set
	607	!-- ug and vg to zero there
	608	ug(nzt+1) = 0.0
	609	vg(nzt+1) = 0.0
	610
	611	ENDIF
	612
	613	!
[1]	614	!-- Calculate virtual potential temperature
[75]	615	IF ( humidity ) vpt = pt * ( 1.0 + 0.61 * q )
[1]	616
	617	!
	618	!-- Store initial profiles for output purposes etc.
	619	hom(:,1,5,:) = SPREAD( u(:,nys,nxl), 2, statistic_regions+1 )
	620	hom(:,1,6,:) = SPREAD( v(:,nys,nxl), 2, statistic_regions+1 )
	621	IF ( ibc_uv_b == 0 ) THEN
	622	hom(nzb,1,5,:) = -hom(nzb+1,1,5,:) ! due to satisfying the Dirichlet
	623	hom(nzb,1,6,:) = -hom(nzb+1,1,6,:) ! condition with an extra layer
	624	! below the surface where the u and v component change their sign
	625	ENDIF
	626	hom(:,1,7,:) = SPREAD( pt(:,nys,nxl), 2, statistic_regions+1 )
	627	hom(:,1,23,:) = SPREAD( km(:,nys,nxl), 2, statistic_regions+1 )
	628	hom(:,1,24,:) = SPREAD( kh(:,nys,nxl), 2, statistic_regions+1 )
	629
[97]	630	IF ( ocean ) THEN
	631	!
	632	!-- Store initial salinity profile
	633	hom(:,1,26,:) = SPREAD( sa(:,nys,nxl), 2, statistic_regions+1 )
	634	ENDIF
[1]	635
[75]	636	IF ( humidity ) THEN
[1]	637	!
	638	!-- Store initial profile of total water content, virtual potential
	639	!-- temperature
	640	hom(:,1,26,:) = SPREAD( q(:,nys,nxl), 2, statistic_regions+1 )
	641	hom(:,1,29,:) = SPREAD( vpt(:,nys,nxl), 2, statistic_regions+1 )
	642	IF ( cloud_physics .OR. cloud_droplets ) THEN
	643	!
	644	!-- Store initial profile of specific humidity and potential
	645	!-- temperature
	646	hom(:,1,27,:) = SPREAD( q(:,nys,nxl), 2, statistic_regions+1 )
	647	hom(:,1,28,:) = SPREAD( pt(:,nys,nxl), 2, statistic_regions+1 )
	648	ENDIF
	649	ENDIF
	650
	651	IF ( passive_scalar ) THEN
	652	!
	653	!-- Store initial scalar profile
	654	hom(:,1,26,:) = SPREAD( q(:,nys,nxl), 2, statistic_regions+1 )
	655	ENDIF
	656
	657	!
[19]	658	!-- Initialize fluxes at bottom surface
[1]	659	IF ( use_surface_fluxes ) THEN
	660
	661	IF ( constant_heatflux ) THEN
	662	!
	663	!-- Heat flux is prescribed
	664	IF ( random_heatflux ) THEN
	665	CALL disturb_heatflux
	666	ELSE
	667	shf = surface_heatflux
	668	!
	669	!-- Over topography surface_heatflux is replaced by wall_heatflux(0)
	670	IF ( TRIM( topography ) /= 'flat' ) THEN
	671	DO i = nxl-1, nxr+1
	672	DO j = nys-1, nyn+1
	673	IF ( nzb_s_inner(j,i) /= 0 ) THEN
	674	shf(j,i) = wall_heatflux(0)
	675	ENDIF
	676	ENDDO
	677	ENDDO
	678	ENDIF
	679	ENDIF
	680	IF ( ASSOCIATED( shf_m ) ) shf_m = shf
	681	ENDIF
	682
	683	!
	684	!-- Determine the near-surface water flux
[75]	685	IF ( humidity .OR. passive_scalar ) THEN
[1]	686	IF ( constant_waterflux ) THEN
	687	qsws = surface_waterflux
	688	IF ( ASSOCIATED( qsws_m ) ) qsws_m = qsws
	689	ENDIF
	690	ENDIF
	691
	692	ENDIF
	693
	694	!
[19]	695	!-- Initialize fluxes at top surface
[94]	696	!-- Currently, only the heatflux and salinity flux can be prescribed.
	697	!-- The latent flux is zero in this case!
[19]	698	IF ( use_top_fluxes ) THEN
	699
	700	IF ( constant_top_heatflux ) THEN
	701	!
	702	!-- Heat flux is prescribed
	703	tswst = top_heatflux
	704	IF ( ASSOCIATED( tswst_m ) ) tswst_m = tswst
	705
[75]	706	IF ( humidity .OR. passive_scalar ) THEN
[19]	707	qswst = 0.0
	708	IF ( ASSOCIATED( qswst_m ) ) qswst_m = qswst
	709	ENDIF
[94]	710
	711	IF ( ocean ) THEN
[95]	712	saswsb = bottom_salinityflux
[94]	713	saswst = top_salinityflux
	714	ENDIF
[102]	715	ENDIF
[19]	716
[102]	717	!
	718	!-- Initialization in case of a coupled model run
	719	IF ( coupling_mode == 'ocean_to_atmosphere' ) THEN
	720	tswst = 0.0
	721	IF ( ASSOCIATED( tswst_m ) ) tswst_m = tswst
	722	ENDIF
	723
[19]	724	ENDIF
	725
	726	!
[1]	727	!-- Initialize Prandtl layer quantities
	728	IF ( prandtl_layer ) THEN
	729
	730	z0 = roughness_length
	731
	732	IF ( .NOT. constant_heatflux ) THEN
	733	!
	734	!-- Surface temperature is prescribed. Here the heat flux cannot be
	735	!-- simply estimated, because therefore rif, u* and theta* would have
	736	!-- to be computed by iteration. This is why the heat flux is assumed
	737	!-- to be zero before the first time step. It approaches its correct
	738	!-- value in the course of the first few time steps.
	739	shf = 0.0
	740	IF ( ASSOCIATED( shf_m ) ) shf_m = 0.0
	741	ENDIF
	742
[75]	743	IF ( humidity .OR. passive_scalar ) THEN
[1]	744	IF ( .NOT. constant_waterflux ) THEN
	745	qsws = 0.0
	746	IF ( ASSOCIATED( qsws_m ) ) qsws_m = 0.0
	747	ENDIF
	748	ENDIF
	749
	750	ENDIF
	751
	752	!
	753	!-- Calculate the initial volume flow at the right and north boundary
	754	IF ( conserve_volume_flow ) THEN
	755
	756	volume_flow_initial_l = 0.0
	757	volume_flow_area_l = 0.0
	758
	759	IF ( nxr == nx ) THEN
	760	DO j = nys, nyn
	761	DO k = nzb_2d(j,nx) + 1, nzt
	762	volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
	763	u(k,j,nx) * dzu(k)
	764	volume_flow_area_l(1) = volume_flow_area_l(1) + dzu(k)
	765	ENDDO
[51]	766	!
	767	!-- Correction if velocity at nzb+1 has been set zero further above
	768	volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
	769	u_nzb_p1_for_vfc(j)
[1]	770	ENDDO
	771	ENDIF
	772
	773	IF ( nyn == ny ) THEN
	774	DO i = nxl, nxr
	775	DO k = nzb_2d(ny,i) + 1, nzt
	776	volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
	777	v(k,ny,i) * dzu(k)
	778	volume_flow_area_l(2) = volume_flow_area_l(2) + dzu(k)
	779	ENDDO
[51]	780	!
	781	!-- Correction if velocity at nzb+1 has been set zero further above
	782	volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
	783	v_nzb_p1_for_vfc(i)
[1]	784	ENDDO
	785	ENDIF
	786
	787	#if defined( __parallel )
	788	CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&
	789	2, MPI_REAL, MPI_SUM, comm2d, ierr )
	790	CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1), &
	791	2, MPI_REAL, MPI_SUM, comm2d, ierr )
	792	#else
	793	volume_flow_initial = volume_flow_initial_l
	794	volume_flow_area = volume_flow_area_l
	795	#endif
	796	ENDIF
	797
	798	!
	799	!-- For the moment, perturbation pressure and vertical velocity are zero
	800	p = 0.0; w = 0.0
	801
	802	!
	803	!-- Initialize array sums (must be defined in first call of pres)
	804	sums = 0.0
	805
	806	!
[72]	807	!-- Treating cloud physics, liquid water content and precipitation amount
	808	!-- are zero at beginning of the simulation
	809	IF ( cloud_physics ) THEN
	810	ql = 0.0
	811	IF ( precipitation ) precipitation_amount = 0.0
	812	ENDIF
[1]	813
	814	!
	815	!-- Initialize spectra
	816	IF ( dt_dosp /= 9999999.9 ) THEN
	817	spectrum_x = 0.0
	818	spectrum_y = 0.0
	819	ENDIF
	820
	821	!
	822	!-- Impose vortex with vertical axis on the initial velocity profile
	823	IF ( INDEX( initializing_actions, 'initialize_vortex' ) /= 0 ) THEN
	824	CALL init_rankine
	825	ENDIF
	826
	827	!
	828	!-- Impose temperature anomaly (advection test only)
	829	IF ( INDEX( initializing_actions, 'initialize_ptanom' ) /= 0 ) THEN
	830	CALL init_pt_anomaly
	831	ENDIF
	832
	833	!
	834	!-- If required, change the surface temperature at the start of the 3D run
	835	IF ( pt_surface_initial_change /= 0.0 ) THEN
	836	pt(nzb,:,:) = pt(nzb,:,:) + pt_surface_initial_change
	837	ENDIF
	838
	839	!
	840	!-- If required, change the surface humidity/scalar at the start of the 3D
	841	!-- run
[75]	842	IF ( ( humidity .OR. passive_scalar ) .AND. &
[1]	843	q_surface_initial_change /= 0.0 ) THEN
	844	q(nzb,:,:) = q(nzb,:,:) + q_surface_initial_change
	845	ENDIF
	846
	847	!
	848	!-- Initialize the random number generator (from numerical recipes)
	849	CALL random_function_ini
	850
	851	!
	852	!-- Impose random perturbation on the horizontal velocity field and then
	853	!-- remove the divergences from the velocity field
	854	IF ( create_disturbances ) THEN
[75]	855	CALL disturb_field( nzb_u_inner, tend, u )
	856	CALL disturb_field( nzb_v_inner, tend, v )
[1]	857	n_sor = nsor_ini
	858	CALL pres
	859	n_sor = nsor
	860	ENDIF
	861
	862	!
	863	!-- Once again set the perturbation pressure explicitly to zero in order to
	864	!-- assure that it does not generate any divergences in the first time step.
	865	!-- At t=0 the velocity field is free of divergence (as constructed above).
	866	!-- Divergences being created during a time step are not yet known and thus
	867	!-- cannot be corrected during the time step yet.
	868	p = 0.0
	869
	870	!
	871	!-- Initialize old and new time levels.
	872	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	873	e_m = e; pt_m = pt; u_m = u; v_m = v; w_m = w; kh_m = kh; km_m = km
	874	ELSE
	875	te_m = 0.0; tpt_m = 0.0; tu_m = 0.0; tv_m = 0.0; tw_m = 0.0
	876	ENDIF
	877	e_p = e; pt_p = pt; u_p = u; v_p = v; w_p = w
	878
[75]	879	IF ( humidity .OR. passive_scalar ) THEN
[1]	880	IF ( ASSOCIATED( q_m ) ) q_m = q
	881	IF ( timestep_scheme(1:5) == 'runge' ) tq_m = 0.0
	882	q_p = q
[75]	883	IF ( humidity .AND. ASSOCIATED( vpt_m ) ) vpt_m = vpt
[1]	884	ENDIF
	885
[94]	886	IF ( ocean ) THEN
	887	tsa_m = 0.0
	888	sa_p = sa
	889	ENDIF
	890
[73]	891	!
	892	!-- Initialize old timelevels needed for radiation boundary conditions
	893	IF ( outflow_l ) THEN
[106]	894	u_m_l(:,:,:) = u(:,:,1:2)
	895	v_m_l(:,:,:) = v(:,:,0:1)
	896	w_m_l(:,:,:) = w(:,:,0:1)
[73]	897	ENDIF
	898	IF ( outflow_r ) THEN
[106]	899	u_m_r(:,:,:) = u(:,:,nx-1:nx)
	900	v_m_r(:,:,:) = v(:,:,nx-1:nx)
	901	w_m_r(:,:,:) = w(:,:,nx-1:nx)
[73]	902	ENDIF
	903	IF ( outflow_s ) THEN
[106]	904	u_m_s(:,:,:) = u(:,0:1,:)
	905	v_m_s(:,:,:) = v(:,1:2,:)
	906	w_m_s(:,:,:) = w(:,0:1,:)
[73]	907	ENDIF
	908	IF ( outflow_n ) THEN
[106]	909	u_m_n(:,:,:) = u(:,ny-1:ny,:)
	910	v_m_n(:,:,:) = v(:,ny-1:ny,:)
	911	w_m_n(:,:,:) = w(:,ny-1:ny,:)
[73]	912	ENDIF
	913
[1]	914	ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data' ) &
	915	THEN
	916	!
	917	!-- Read binary data from restart file
	918	CALL read_3d_binary
	919
	920	!
	921	!-- Calculate initial temperature field and other constants used in case
	922	!-- of a sloping surface
	923	IF ( sloping_surface ) CALL init_slope
	924
	925	!
	926	!-- Initialize new time levels (only done in order to set boundary values
	927	!-- including ghost points)
	928	e_p = e; pt_p = pt; u_p = u; v_p = v; w_p = w
[75]	929	IF ( humidity .OR. passive_scalar ) q_p = q
[94]	930	IF ( ocean ) sa_p = sa
[1]	931
	932	ELSE
	933	!
	934	!-- Actually this part of the programm should not be reached
	935	IF ( myid == 0 ) PRINT*,'+++ init_3d_model: unknown initializing ', &
	936	'problem'
	937	CALL local_stop
	938	ENDIF
	939
	940	!
[138]	941	!-- Initialization of the leaf area density
	942	IF ( plant_canopy ) THEN
	943
	944	SELECT CASE ( TRIM( canopy_mode ) )
	945
	946	CASE( 'block' )
	947
	948	DO i = nxl-1, nxr+1
	949	DO j = nys-1, nyn+1
	950	lad_s(:,j,i) = lad(:)
	951	cdc(:,j,i) = drag_coefficient
	952	ENDDO
	953	ENDDO
	954
	955	CASE DEFAULT
	956
	957	!
	958	!-- The DEFAULT case is reached either if the parameter
	959	!-- canopy mode contains a wrong character string or if the
	960	!-- user has coded a special case in the user interface.
	961	!-- There, the subroutine user_init_plant_canopy checks
	962	!-- which of these two conditions applies.
	963	CALL user_init_plant_canopy
	964
	965	END SELECT
	966
	967	CALL exchange_horiz( lad_s )
	968	CALL exchange_horiz( cdc )
	969
	970	DO i = nxl, nxr
	971	DO j = nys, nyn
	972	DO k = nzb, nzt+1
	973	lad_u(k,j,i) = 0.5 * ( lad_s(k,j,i-1) + lad_s(k,j,i) )
	974	lad_v(k,j,i) = 0.5 * ( lad_s(k,j-1,i) + lad_s(k,j,i) )
	975	ENDDO
	976	DO k = nzb, nzt
	977	lad_w(k,j,i) = 0.5 * ( lad_s(k+1,j,i) + lad_s(k,j,i) )
	978	ENDDO
	979	ENDDO
	980	ENDDO
	981
	982	lad_w(nzt+1,:,:) = lad_w(nzt,:,:)
	983
	984	CALL exchange_horiz( lad_u )
	985	CALL exchange_horiz( lad_v )
	986	CALL exchange_horiz( lad_w )
	987
	988	ENDIF
	989
	990	!
[1]	991	!-- If required, initialize dvrp-software
	992	IF ( dt_dvrp /= 9999999.9 ) CALL init_dvrp
	993
[96]	994	IF ( ocean ) THEN
[1]	995	!
[96]	996	!-- Initialize quantities needed for the ocean model
	997	CALL init_ocean
	998	ELSE
	999	!
	1000	!-- Initialize quantities for handling cloud physics
	1001	!-- This routine must be called before init_particles, because
	1002	!-- otherwise, array pt_d_t, needed in data_output_dvrp (called by
	1003	!-- init_particles) is not defined.
	1004	CALL init_cloud_physics
	1005	ENDIF
[1]	1006
	1007	!
	1008	!-- If required, initialize particles
[63]	1009	IF ( particle_advection ) CALL init_particles
[1]	1010
	1011	!
	1012	!-- Initialize quantities for special advections schemes
	1013	CALL init_advec
	1014
	1015	!
	1016	!-- Initialize Rayleigh damping factors
	1017	rdf = 0.0
	1018	IF ( rayleigh_damping_factor /= 0.0 ) THEN
[108]	1019	IF ( .NOT. ocean ) THEN
	1020	DO k = nzb+1, nzt
	1021	IF ( zu(k) >= rayleigh_damping_height ) THEN
	1022	rdf(k) = rayleigh_damping_factor * &
[1]	1023	( SIN( pi * 0.5 * ( zu(k) - rayleigh_damping_height ) &
	1024	/ ( zu(nzt) - rayleigh_damping_height ) )&
	1025	)**2
[108]	1026	ENDIF
	1027	ENDDO
	1028	ELSE
	1029	DO k = nzt, nzb+1, -1
	1030	IF ( zu(k) <= rayleigh_damping_height ) THEN
	1031	rdf(k) = rayleigh_damping_factor * &
	1032	( SIN( pi * 0.5 * ( rayleigh_damping_height - zu(k) ) &
	1033	/ ( rayleigh_damping_height - zu(nzb+1)))&
	1034	)**2
	1035	ENDIF
	1036	ENDDO
	1037	ENDIF
[1]	1038	ENDIF
	1039
	1040	!
	1041	!-- Initialize diffusivities used within the outflow damping layer in case of
	1042	!-- non-cyclic lateral boundaries. A linear increase is assumed over the first
	1043	!-- half of the width of the damping layer
[73]	1044	IF ( bc_lr == 'dirichlet/radiation' ) THEN
[1]	1045
	1046	DO i = nxl-1, nxr+1
[73]	1047	IF ( i >= nx - outflow_damping_width ) THEN
	1048	km_damp_x(i) = km_damp_max * MIN( 1.0, &
	1049	( i - ( nx - outflow_damping_width ) ) / &
	1050	REAL( outflow_damping_width/2 ) &
	1051	)
	1052	ELSE
	1053	km_damp_x(i) = 0.0
	1054	ENDIF
	1055	ENDDO
[1]	1056
[73]	1057	ELSEIF ( bc_lr == 'radiation/dirichlet' ) THEN
[1]	1058
[73]	1059	DO i = nxl-1, nxr+1
	1060	IF ( i <= outflow_damping_width ) THEN
	1061	km_damp_x(i) = km_damp_max * MIN( 1.0, &
	1062	( outflow_damping_width - i ) / &
	1063	REAL( outflow_damping_width/2 ) &
	1064	)
	1065	ELSE
	1066	km_damp_x(i) = 0.0
	1067	ENDIF
	1068	ENDDO
[1]	1069
[73]	1070	ENDIF
[1]	1071
[73]	1072	IF ( bc_ns == 'dirichlet/radiation' ) THEN
[1]	1073
[73]	1074	DO j = nys-1, nyn+1
	1075	IF ( j >= ny - outflow_damping_width ) THEN
	1076	km_damp_y(j) = km_damp_max * MIN( 1.0, &
	1077	( j - ( ny - outflow_damping_width ) ) / &
	1078	REAL( outflow_damping_width/2 ) &
	1079	)
	1080	ELSE
	1081	km_damp_y(j) = 0.0
[1]	1082	ENDIF
	1083	ENDDO
	1084
[73]	1085	ELSEIF ( bc_ns == 'radiation/dirichlet' ) THEN
[1]	1086
	1087	DO j = nys-1, nyn+1
[73]	1088	IF ( j <= outflow_damping_width ) THEN
	1089	km_damp_y(j) = km_damp_max * MIN( 1.0, &
	1090	( outflow_damping_width - j ) / &
	1091	REAL( outflow_damping_width/2 ) &
	1092	)
	1093	ELSE
	1094	km_damp_y(j) = 0.0
[1]	1095	ENDIF
[73]	1096	ENDDO
[1]	1097
	1098	ENDIF
	1099
	1100	!
	1101	!-- Initialize local summation arrays for UP flow_statistics. This is necessary
	1102	!-- because they may not yet have been initialized when they are called from
	1103	!-- flow_statistics (or - depending on the chosen model run - are never
	1104	!-- initialized)
	1105	sums_divnew_l = 0.0
	1106	sums_divold_l = 0.0
	1107	sums_l_l = 0.0
	1108	sums_up_fraction_l = 0.0
	1109	sums_wsts_bc_l = 0.0
	1110
	1111	!
	1112	!-- Pre-set masks for regional statistics. Default is the total model domain.
	1113	rmask = 1.0
	1114
	1115	!
[51]	1116	!-- User-defined initializing actions. Check afterwards, if maximum number
	1117	!-- of allowed timeseries is not exceeded
[1]	1118	CALL user_init
	1119
[51]	1120	IF ( dots_num > dots_max ) THEN
	1121	IF ( myid == 0 ) THEN
	1122	PRINT*, '+++ user_init: number of time series quantities exceeds', &
	1123	' its maximum of dots_max = ', dots_max
	1124	PRINT*, ' Please increase dots_max in modules.f90.'
	1125	ENDIF
	1126	CALL local_stop
	1127	ENDIF
	1128
[1]	1129	!
	1130	!-- Input binary data file is not needed anymore. This line must be placed
	1131	!-- after call of user_init!
	1132	CALL close_file( 13 )
	1133
	1134	!
	1135	!-- Compute total sum of active mask grid points
	1136	!-- ngp_2dh: number of grid points of a horizontal cross section through the
	1137	!-- total domain
	1138	!-- ngp_3d: number of grid points of the total domain
[132]	1139	ngp_2dh_outer_l = 0
	1140	ngp_2dh_outer = 0
	1141	ngp_2dh_s_inner_l = 0
	1142	ngp_2dh_s_inner = 0
	1143	ngp_2dh_l = 0
	1144	ngp_2dh = 0
	1145	ngp_3d_inner_l = 0
	1146	ngp_3d_inner = 0
	1147	ngp_3d = 0
	1148	ngp_sums = ( nz + 2 ) * ( pr_palm + max_pr_user )
[1]	1149
	1150	DO sr = 0, statistic_regions
	1151	DO i = nxl, nxr
	1152	DO j = nys, nyn
	1153	IF ( rmask(j,i,sr) == 1.0 ) THEN
	1154	!
	1155	!-- All xy-grid points
	1156	ngp_2dh_l(sr) = ngp_2dh_l(sr) + 1
	1157	!
	1158	!-- xy-grid points above topography
	1159	DO k = nzb_s_outer(j,i), nz + 1
	1160	ngp_2dh_outer_l(k,sr) = ngp_2dh_outer_l(k,sr) + 1
	1161	ENDDO
[132]	1162	DO k = nzb_s_inner(j,i), nz + 1
	1163	ngp_2dh_s_inner_l(k,sr) = ngp_2dh_s_inner_l(k,sr) + 1
	1164	ENDDO
[1]	1165	!
	1166	!-- All grid points of the total domain above topography
	1167	ngp_3d_inner_l(sr) = ngp_3d_inner_l(sr) + &
	1168	( nz - nzb_s_inner(j,i) + 2 )
	1169	ENDIF
	1170	ENDDO
	1171	ENDDO
	1172	ENDDO
	1173
	1174	sr = statistic_regions + 1
	1175	#if defined( __parallel )
	1176	CALL MPI_ALLREDUCE( ngp_2dh_l(0), ngp_2dh(0), sr, MPI_INTEGER, MPI_SUM, &
	1177	comm2d, ierr )
	1178	CALL MPI_ALLREDUCE( ngp_2dh_outer_l(0,0), ngp_2dh_outer(0,0), (nz+2)*sr, &
	1179	MPI_INTEGER, MPI_SUM, comm2d, ierr )
[132]	1180	CALL MPI_ALLREDUCE( ngp_2dh_s_inner_l(0,0), ngp_2dh_s_inner(0,0), &
	1181	(nz+2)*sr, MPI_INTEGER, MPI_SUM, comm2d, ierr )
[1]	1182	CALL MPI_ALLREDUCE( ngp_3d_inner_l(0), ngp_3d_inner(0), sr, MPI_INTEGER, &
	1183	MPI_SUM, comm2d, ierr )
	1184	#else
[132]	1185	ngp_2dh = ngp_2dh_l
	1186	ngp_2dh_outer = ngp_2dh_outer_l
	1187	ngp_2dh_s_inner = ngp_2dh_s_inner_l
	1188	ngp_3d_inner = ngp_3d_inner_l
[1]	1189	#endif
	1190
	1191	ngp_3d = ngp_2dh * ( nz + 2 )
	1192
	1193	!
	1194	!-- Set a lower limit of 1 in order to avoid zero divisions in flow_statistics,
	1195	!-- buoyancy, etc. A zero value will occur for cases where all grid points of
	1196	!-- the respective subdomain lie below the surface topography
	1197	ngp_2dh_outer = MAX( 1, ngp_2dh_outer(:,:) )
	1198	ngp_3d_inner = MAX( 1, ngp_3d_inner(:) )
	1199
	1200	DEALLOCATE( ngp_2dh_l, ngp_2dh_outer_l, ngp_3d_inner_l )
	1201
	1202
	1203	END SUBROUTINE init_3d_model

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

Context Navigation

source: palm/trunk/SOURCE/init_3d_model.f90 @ 141

Download in other formats: