Home

Context Navigation

init_3d_model.f90 @ 707

Last change on this file since 707 was 707, checked in by raasch, 13 years ago

New:
---

In case of multigrid method, on coarse grid levels, gathered data are
identically processed on all PEs (before, on PE0 only), so that the subsequent
scattering of data is not neccessary any more. (modules, init_pegrid, poismg)

Changed:

Calculation of weighted average of p is now handled in the same way
regardless of the number of ghost layers (advection scheme). (pres)

multigrid and sor method are using p_loc for iterative
advancements of pressure. p_sub removed. (init_3d_model, modules, poismg, pres, sor)

bc_lr and bc_ns replaced by bc_lr_dirrad, bc_lr_raddir, bc_ns_dirrad, bc_ns_raddir
for speed optimization. (calc_spectra, check_parameters, exchange_horiz,
exchange_horiz_2d, header, init_3d_model, init_grid, init_pegrid, modules,
poismg, pres, sor, time_integration, timestep)

grid_level directly used as index for MPI data type arrays. (exchange_horiz,
poismg)

initial assignments of zero to array p for iterative solvers only (init_3d_model)

Errors:

localsum calculation modified for proper OpenMP reduction. (pres)

Bugfix: bottom (nzb) and top (nzt+1) boundary conditions set in routines
resid and restrict. They were missed before, which may have led to
unpredictable results. (poismg)

Property svn:keywords set to Id

File size: 55.3 KB

Rev	Line
[1]	1	#if defined( __ibmy_special )
	2	@PROCESS NOOPTimize
	3	#endif
	4	SUBROUTINE init_3d_model
	5
	6	!------------------------------------------------------------------------------!
[254]	7	! Current revisions:
[1]	8	! -----------------
[707]	9	! p_sub renamed p_loc and allocated depending on the chosen pressure solver,
	10	! initial assignments of zero to array p for iterative solvers only,
	11	! bc_lr/ns replaced by bc_lr/ns_dirrad/raddir
[674]	12	!
[668]	13	! Former revisions:
	14	! -----------------
[681]	15	! 680 2011-02-04 23:16:06Z gryschka $
	16	! bugfix: volume_flow_control
[668]	17	!
[674]	18	! 673 2011-01-18 16:19:48Z suehring
	19	! weight_substep (moved from advec_ws) and weight_pres added.
	20	! Allocate p_sub when using Multigrid or SOR solver.
	21	! Call of ws_init moved behind the if requests.
	22	!
[668]	23	! 667 2010-12-23 12:06:00Z suehring/gryschka
[667]	24	! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
	25	! allocation of arrays. Calls of exchange_horiz are modified.
	26	! Call ws_init to initialize arrays needed for statistical evaluation and
	27	! optimization when ws-scheme is used.
	28	! Initial volume flow is now calculated by using the variable hom_sum.
	29	! Therefore the correction of initial volume flow for non-flat topography
	30	! removed (removed u_nzb_p1_for_vfc and v_nzb_p1_for_vfc)
	31	! Changed surface boundary conditions for u and v in case of ibc_uv_b == 0 from
	32	! mirror bc to dirichlet boundary conditions (u=v=0), so that k=nzb is
	33	! representative for the height z0
	34	! Bugfix: type conversion of '1' to 64bit for the MAX function (ngp_3d_inner)
	35	!
[623]	36	! 622 2010-12-10 08:08:13Z raasch
	37	! optional barriers included in order to speed up collective operations
	38	!
[561]	39	! 560 2010-09-09 10:06:09Z weinreis
	40	! bugfix: correction of calculating ngp_3d for 64 bit
	41	!
[486]	42	! 485 2010-02-05 10:57:51Z raasch
	43	! calculation of ngp_3d + ngp_3d_inner changed because they have now 64 bit
	44	!
[482]	45	! 407 2009-12-01 15:01:15Z maronga
	46	! var_ts is replaced by dots_max
	47	! Enabled passive scalar/humidity wall fluxes for non-flat topography
	48	!
[392]	49	! 388 2009-09-23 09:40:33Z raasch
[388]	50	! Initialization of prho added.
[359]	51	! bugfix: correction of initial volume flow for non-flat topography
	52	! bugfix: zero initialization of arrays within buildings for 'cyclic_fill'
[333]	53	! bugfix: avoid that ngp_2dh_s_inner becomes zero
[328]	54	! initializing_actions='read_data_for_recycling' renamed to 'cyclic_fill', now
	55	! independent of turbulent_inflow
[254]	56	! Output of messages replaced by message handling routine.
[240]	57	! Set the starting level and the vertical smoothing factor used for
	58	! the external pressure gradient
[254]	59	! +conserve_volume_flow_mode: 'default', 'initial_profiles', 'inflow_profile'
[241]	60	! and 'bulk_velocity'
[292]	61	! If the inversion height calculated by the prerun is zero,
	62	! inflow_damping_height must be explicitly specified.
[139]	63	!
[198]	64	! 181 2008-07-30 07:07:47Z raasch
	65	! bugfix: zero assignments to tendency arrays in case of restarts,
	66	! further extensions and modifications in the initialisation of the plant
	67	! canopy model,
	68	! allocation of hom_sum moved to parin, initialization of spectrum_x\|y directly
	69	! after allocating theses arrays,
	70	! read data for recycling added as new initialization option,
	71	! dummy allocation for diss
	72	!
[139]	73	! 138 2007-11-28 10:03:58Z letzel
[132]	74	! New counter ngp_2dh_s_inner.
	75	! Allow new case bc_uv_t = 'dirichlet_0' for channel flow.
	76	! Corrected calculation of initial volume flow for 'set_1d-model_profiles' and
	77	! 'set_constant_profiles' in case of buildings in the reference cross-sections.
[77]	78	!
[110]	79	! 108 2007-08-24 15:10:38Z letzel
	80	! Flux initialization in case of coupled runs, +momentum fluxes at top boundary,
	81	! +arrays for phase speed c_u, c_v, c_w, indices for u\|v\|w_m_l\|r changed
	82	! +qswst_remote in case of atmosphere model with humidity coupled to ocean
	83	! Rayleigh damping for ocean, optionally calculate km and kh from initial
	84	! TKE e_init
	85	!
[98]	86	! 97 2007-06-21 08:23:15Z raasch
	87	! Initialization of salinity, call of init_ocean
	88	!
[90]	89	! 87 2007-05-22 15:46:47Z raasch
	90	! var_hom and var_sum renamed pr_palm
	91	!
[77]	92	! 75 2007-03-22 09:54:05Z raasch
[73]	93	! Arrays for radiation boundary conditions are allocated (u_m_l, u_m_r, etc.),
	94	! bugfix for cases with the outflow damping layer extending over more than one
[75]	95	! subdomain, moisture renamed humidity,
	96	! new initializing action "by_user" calls user_init_3d_model,
[72]	97	! precipitation_amount/rate, ts_value are allocated, +module netcdf_control,
[51]	98	! initial velocities at nzb+1 are regarded for volume
	99	! flow control in case they have been set zero before (to avoid small timesteps)
[75]	100	! -uvmean_outflow, uxrp, vynp eliminated
[1]	101	!
[39]	102	! 19 2007-02-23 04:53:48Z raasch
	103	! +handling of top fluxes
	104	!
[3]	105	! RCS Log replace by Id keyword, revision history cleaned up
	106	!
[1]	107	! Revision 1.49 2006/08/22 15:59:07 raasch
	108	! No optimization of this file on the ibmy (Yonsei Univ.)
	109	!
	110	! Revision 1.1 1998/03/09 16:22:22 raasch
	111	! Initial revision
	112	!
	113	!
	114	! Description:
	115	! ------------
	116	! Allocation of arrays and initialization of the 3D model via
	117	! a) pre-run the 1D model
	118	! or
	119	! b) pre-set constant linear profiles
	120	! or
	121	! c) read values of a previous run
	122	!------------------------------------------------------------------------------!
	123
[667]	124	USE advec_ws
[1]	125	USE arrays_3d
	126	USE averaging
[72]	127	USE cloud_parameters
[1]	128	USE constants
	129	USE control_parameters
	130	USE cpulog
	131	USE indices
	132	USE interfaces
	133	USE model_1d
[51]	134	USE netcdf_control
[1]	135	USE particle_attributes
	136	USE pegrid
	137	USE profil_parameter
	138	USE random_function_mod
	139	USE statistics
	140
	141	IMPLICIT NONE
	142
[559]	143	INTEGER :: i, ind_array(1), j, k, sr
[1]	144
[485]	145	INTEGER, DIMENSION(:), ALLOCATABLE :: ngp_2dh_l
[1]	146
[132]	147	INTEGER, DIMENSION(:,:), ALLOCATABLE :: ngp_2dh_outer_l, &
	148	ngp_2dh_s_inner_l
[1]	149
[153]	150	REAL :: a, b
	151
[1]	152	REAL, DIMENSION(1:2) :: volume_flow_area_l, volume_flow_initial_l
	153
[485]	154	REAL, DIMENSION(:), ALLOCATABLE :: ngp_3d_inner_l, ngp_3d_inner_tmp
[1]	155
[485]	156
[1]	157	!
	158	!-- Allocate arrays
	159	ALLOCATE( ngp_2dh(0:statistic_regions), ngp_2dh_l(0:statistic_regions), &
	160	ngp_3d(0:statistic_regions), &
	161	ngp_3d_inner(0:statistic_regions), &
	162	ngp_3d_inner_l(0:statistic_regions), &
[485]	163	ngp_3d_inner_tmp(0:statistic_regions), &
[1]	164	sums_divnew_l(0:statistic_regions), &
	165	sums_divold_l(0:statistic_regions) )
[240]	166	ALLOCATE( dp_smooth_factor(nzb:nzt), rdf(nzb+1:nzt) )
[143]	167	ALLOCATE( ngp_2dh_outer(nzb:nzt+1,0:statistic_regions), &
[1]	168	ngp_2dh_outer_l(nzb:nzt+1,0:statistic_regions), &
[132]	169	ngp_2dh_s_inner(nzb:nzt+1,0:statistic_regions), &
	170	ngp_2dh_s_inner_l(nzb:nzt+1,0:statistic_regions), &
[667]	171	rmask(nysg:nyng,nxlg:nxrg,0:statistic_regions), &
[87]	172	sums(nzb:nzt+1,pr_palm+max_pr_user), &
	173	sums_l(nzb:nzt+1,pr_palm+max_pr_user,0:threads_per_task-1), &
[1]	174	sums_l_l(nzb:nzt+1,0:statistic_regions,0:threads_per_task-1), &
	175	sums_up_fraction_l(10,3,0:statistic_regions), &
[48]	176	sums_wsts_bc_l(nzb:nzt+1,0:statistic_regions), &
[394]	177	ts_value(dots_max,0:statistic_regions) )
[667]	178	ALLOCATE( km_damp_x(nxlg:nxrg), km_damp_y(nysg:nyng) )
[1]	179
[667]	180	ALLOCATE( rif_1(nysg:nyng,nxlg:nxrg), shf_1(nysg:nyng,nxlg:nxrg), &
	181	ts(nysg:nyng,nxlg:nxrg), tswst_1(nysg:nyng,nxlg:nxrg), &
	182	us(nysg:nyng,nxlg:nxrg), usws_1(nysg:nyng,nxlg:nxrg), &
	183	uswst_1(nysg:nyng,nxlg:nxrg), &
	184	vsws_1(nysg:nyng,nxlg:nxrg), &
	185	vswst_1(nysg:nyng,nxlg:nxrg), z0(nysg:nyng,nxlg:nxrg) )
[1]	186
	187	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	188	!
	189	!-- Leapfrog scheme needs two timelevels of diffusion quantities
[667]	190	ALLOCATE( rif_2(nysg:nyng,nxlg:nxrg), &
	191	shf_2(nysg:nyng,nxlg:nxrg), &
	192	tswst_2(nysg:nyng,nxlg:nxrg), &
	193	usws_2(nysg:nyng,nxlg:nxrg), &
	194	uswst_2(nysg:nyng,nxlg:nxrg), &
	195	vswst_2(nysg:nyng,nxlg:nxrg), &
	196	vsws_2(nysg:nyng,nxlg:nxrg) )
[1]	197	ENDIF
	198
[75]	199	ALLOCATE( d(nzb+1:nzta,nys:nyna,nxl:nxra), &
[667]	200	e_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	201	e_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	202	e_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	203	kh_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	204	km_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	205	p(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	206	pt_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	207	pt_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	208	pt_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	209	tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	210	u_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	211	u_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	212	u_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	213	v_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	214	v_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	215	v_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	216	w_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	217	w_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	218	w_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[673]	219	!
[707]	220	!-- Following array is required for perturbation pressure within the iterative
	221	!-- pressure solvers. For the multistep schemes (Runge-Kutta), array p holds
	222	!-- the weighted average of the substeps and cannot be used in the Poisson
	223	!-- solver.
	224	IF ( psolver == 'sor' ) THEN
	225	ALLOCATE( p_loc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	226	ELSEIF ( psolver == 'multigrid' ) THEN
	227	!
	228	!-- For performance reasons, multigrid is using one ghost layer only
	229	ALLOCATE( p_loc(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
[673]	230	ENDIF
[1]	231
	232	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
[667]	233	ALLOCATE( kh_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	234	km_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]	235	ENDIF
	236
[75]	237	IF ( humidity .OR. passive_scalar ) THEN
[1]	238	!
[75]	239	!-- 2D-humidity/scalar arrays
[667]	240	ALLOCATE ( qs(nysg:nyng,nxlg:nxrg), &
	241	qsws_1(nysg:nyng,nxlg:nxrg), &
	242	qswst_1(nysg:nyng,nxlg:nxrg) )
[1]	243
	244	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
[667]	245	ALLOCATE( qsws_2(nysg:nyng,nxlg:nxrg), &
	246	qswst_2(nysg:nyng,nxlg:nxrg) )
[1]	247	ENDIF
	248	!
[75]	249	!-- 3D-humidity/scalar arrays
[667]	250	ALLOCATE( q_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	251	q_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	252	q_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]	253
	254	!
[75]	255	!-- 3D-arrays needed for humidity only
	256	IF ( humidity ) THEN
[667]	257	ALLOCATE( vpt_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]	258
	259	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
[667]	260	ALLOCATE( vpt_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]	261	ENDIF
	262
	263	IF ( cloud_physics ) THEN
	264	!
	265	!-- Liquid water content
[667]	266	ALLOCATE ( ql_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[72]	267	!
	268	!-- Precipitation amount and rate (only needed if output is switched)
[667]	269	ALLOCATE( precipitation_amount(nysg:nyng,nxlg:nxrg), &
	270	precipitation_rate(nysg:nyng,nxlg:nxrg) )
[1]	271	ENDIF
	272
	273	IF ( cloud_droplets ) THEN
	274	!
	275	!-- Liquid water content, change in liquid water content,
	276	!-- real volume of particles (with weighting), volume of particles
[667]	277	ALLOCATE ( ql_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	278	ql_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	279	ql_v(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	280	ql_vp(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]	281	ENDIF
	282
	283	ENDIF
	284
	285	ENDIF
	286
[94]	287	IF ( ocean ) THEN
[667]	288	ALLOCATE( saswsb_1(nysg:nyng,nxlg:nxrg), &
	289	saswst_1(nysg:nyng,nxlg:nxrg) )
	290	ALLOCATE( prho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	291	rho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	292	sa_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	293	sa_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	294	sa_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[388]	295	prho => prho_1
	296	rho => rho_1 ! routines calc_mean_profile and diffusion_e require
	297	! density to be apointer
[108]	298	IF ( humidity_remote ) THEN
[667]	299	ALLOCATE( qswst_remote(nysg:nyng,nxlg:nxrg))
[108]	300	qswst_remote = 0.0
	301	ENDIF
[94]	302	ENDIF
	303
[1]	304	!
	305	!-- 3D-array for storing the dissipation, needed for calculating the sgs
	306	!-- particle velocities
	307	IF ( use_sgs_for_particles ) THEN
[667]	308	ALLOCATE ( diss(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[181]	309	ELSE
	310	ALLOCATE ( diss(2,2,2) ) ! required because diss is used as a
	311	! formal parameter
[1]	312	ENDIF
	313
	314	IF ( dt_dosp /= 9999999.9 ) THEN
	315	ALLOCATE( spectrum_x( 1:nx/2, 1:10, 1:10 ), &
	316	spectrum_y( 1:ny/2, 1:10, 1:10 ) )
[146]	317	spectrum_x = 0.0
	318	spectrum_y = 0.0
[1]	319	ENDIF
	320
	321	!
[138]	322	!-- 3D-arrays for the leaf area density and the canopy drag coefficient
	323	IF ( plant_canopy ) THEN
[667]	324	ALLOCATE ( lad_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	325	lad_u(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	326	lad_v(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	327	lad_w(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	328	cdc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[153]	329
	330	IF ( passive_scalar ) THEN
[667]	331	ALLOCATE ( sls(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	332	sec(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[153]	333	ENDIF
	334
	335	IF ( cthf /= 0.0 ) THEN
[667]	336	ALLOCATE ( lai(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	337	canopy_heat_flux(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[153]	338	ENDIF
	339
[138]	340	ENDIF
	341
	342	!
[51]	343	!-- 4D-array for storing the Rif-values at vertical walls
	344	IF ( topography /= 'flat' ) THEN
[667]	345	ALLOCATE( rif_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg,1:4) )
[51]	346	rif_wall = 0.0
	347	ENDIF
	348
	349	!
[106]	350	!-- Arrays to store velocity data from t-dt and the phase speeds which
	351	!-- are needed for radiation boundary conditions
[73]	352	IF ( outflow_l ) THEN
[667]	353	ALLOCATE( u_m_l(nzb:nzt+1,nysg:nyng,1:2), &
	354	v_m_l(nzb:nzt+1,nysg:nyng,0:1), &
	355	w_m_l(nzb:nzt+1,nysg:nyng,0:1) )
[73]	356	ENDIF
	357	IF ( outflow_r ) THEN
[667]	358	ALLOCATE( u_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx), &
	359	v_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx), &
	360	w_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx) )
[73]	361	ENDIF
[106]	362	IF ( outflow_l .OR. outflow_r ) THEN
[667]	363	ALLOCATE( c_u(nzb:nzt+1,nysg:nyng), c_v(nzb:nzt+1,nysg:nyng), &
	364	c_w(nzb:nzt+1,nysg:nyng) )
[106]	365	ENDIF
[73]	366	IF ( outflow_s ) THEN
[667]	367	ALLOCATE( u_m_s(nzb:nzt+1,0:1,nxlg:nxrg), &
	368	v_m_s(nzb:nzt+1,1:2,nxlg:nxrg), &
	369	w_m_s(nzb:nzt+1,0:1,nxlg:nxrg) )
[73]	370	ENDIF
	371	IF ( outflow_n ) THEN
[667]	372	ALLOCATE( u_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg), &
	373	v_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg), &
	374	w_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg) )
[73]	375	ENDIF
[106]	376	IF ( outflow_s .OR. outflow_n ) THEN
[667]	377	ALLOCATE( c_u(nzb:nzt+1,nxlg:nxrg), c_v(nzb:nzt+1,nxlg:nxrg), &
	378	c_w(nzb:nzt+1,nxlg:nxrg) )
[106]	379	ENDIF
[73]	380
	381	!
[1]	382	!-- Initial assignment of the pointers
	383	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	384
[19]	385	rif_m => rif_1; rif => rif_2
	386	shf_m => shf_1; shf => shf_2
	387	tswst_m => tswst_1; tswst => tswst_2
	388	usws_m => usws_1; usws => usws_2
[102]	389	uswst_m => uswst_1; uswst => uswst_2
[19]	390	vsws_m => vsws_1; vsws => vsws_2
[102]	391	vswst_m => vswst_1; vswst => vswst_2
[1]	392	e_m => e_1; e => e_2; e_p => e_3; te_m => e_3
	393	kh_m => kh_1; kh => kh_2
	394	km_m => km_1; km => km_2
	395	pt_m => pt_1; pt => pt_2; pt_p => pt_3; tpt_m => pt_3
	396	u_m => u_1; u => u_2; u_p => u_3; tu_m => u_3
	397	v_m => v_1; v => v_2; v_p => v_3; tv_m => v_3
	398	w_m => w_1; w => w_2; w_p => w_3; tw_m => w_3
	399
[75]	400	IF ( humidity .OR. passive_scalar ) THEN
[19]	401	qsws_m => qsws_1; qsws => qsws_2
	402	qswst_m => qswst_1; qswst => qswst_2
[1]	403	q_m => q_1; q => q_2; q_p => q_3; tq_m => q_3
[75]	404	IF ( humidity ) vpt_m => vpt_1; vpt => vpt_2
[1]	405	IF ( cloud_physics ) ql => ql_1
	406	IF ( cloud_droplets ) THEN
	407	ql => ql_1
	408	ql_c => ql_2
	409	ENDIF
	410	ENDIF
	411
	412	ELSE
	413
[19]	414	rif => rif_1
	415	shf => shf_1
	416	tswst => tswst_1
	417	usws => usws_1
[102]	418	uswst => uswst_1
[19]	419	vsws => vsws_1
[102]	420	vswst => vswst_1
[19]	421	e => e_1; e_p => e_2; te_m => e_3; e_m => e_3
	422	kh => kh_1
	423	km => km_1
	424	pt => pt_1; pt_p => pt_2; tpt_m => pt_3; pt_m => pt_3
	425	u => u_1; u_p => u_2; tu_m => u_3; u_m => u_3
	426	v => v_1; v_p => v_2; tv_m => v_3; v_m => v_3
	427	w => w_1; w_p => w_2; tw_m => w_3; w_m => w_3
[1]	428
[75]	429	IF ( humidity .OR. passive_scalar ) THEN
[1]	430	qsws => qsws_1
[19]	431	qswst => qswst_1
[94]	432	q => q_1; q_p => q_2; tq_m => q_3; q_m => q_3
[75]	433	IF ( humidity ) vpt => vpt_1
[1]	434	IF ( cloud_physics ) ql => ql_1
	435	IF ( cloud_droplets ) THEN
	436	ql => ql_1
	437	ql_c => ql_2
	438	ENDIF
	439	ENDIF
	440
[94]	441	IF ( ocean ) THEN
[95]	442	saswsb => saswsb_1
[94]	443	saswst => saswst_1
	444	sa => sa_1; sa_p => sa_2; tsa_m => sa_3
	445	ENDIF
	446
[1]	447	ENDIF
[673]	448
[1]	449	!
[673]	450	!-- Allocate arrays containing the RK coefficient for right evaluation of
	451	!-- perturbation pressure and turbulent fluxes. At this point it is needed
	452	!-- for right pressure correction during initialization. Further below
	453	!-- the real values will be set.
	454	ALLOCATE (weight_substep(1:intermediate_timestep_count_max), &
	455	weight_pres(1:intermediate_timestep_count_max) )
	456	weight_substep = 1.
	457	weight_pres = 1.
	458	intermediate_timestep_count = 1 ! needed for simulated_time=0
	459
	460	!
[1]	461	!-- Initialize model variables
[147]	462	IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. &
[328]	463	TRIM( initializing_actions ) /= 'cyclic_fill' ) THEN
[1]	464	!
	465	!-- First model run of a possible job queue.
	466	!-- Initial profiles of the variables must be computes.
	467	IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 ) THEN
	468	!
	469	!-- Use solutions of the 1D model as initial profiles,
	470	!-- start 1D model
	471	CALL init_1d_model
	472	!
	473	!-- Transfer initial profiles to the arrays of the 3D model
[667]	474	DO i = nxlg, nxrg
	475	DO j = nysg, nyng
[1]	476	e(:,j,i) = e1d
	477	kh(:,j,i) = kh1d
	478	km(:,j,i) = km1d
	479	pt(:,j,i) = pt_init
	480	u(:,j,i) = u1d
	481	v(:,j,i) = v1d
	482	ENDDO
	483	ENDDO
	484
[75]	485	IF ( humidity .OR. passive_scalar ) THEN
[667]	486	DO i = nxlg, nxrg
	487	DO j = nysg, nyng
[1]	488	q(:,j,i) = q_init
	489	ENDDO
	490	ENDDO
	491	ENDIF
	492
	493	IF ( .NOT. constant_diffusion ) THEN
[667]	494	DO i = nxlg, nxrg
	495	DO j = nysg, nyng
[1]	496	e(:,j,i) = e1d
	497	ENDDO
	498	ENDDO
	499	!
	500	!-- Store initial profiles for output purposes etc.
	501	hom(:,1,25,:) = SPREAD( l1d, 2, statistic_regions+1 )
	502
	503	IF ( prandtl_layer ) THEN
	504	rif = rif1d(nzb+1)
	505	ts = 0.0 ! could actually be computed more accurately in the
	506	! 1D model. Update when opportunity arises.
	507	us = us1d
	508	usws = usws1d
	509	vsws = vsws1d
	510	ELSE
	511	ts = 0.0 ! must be set, because used in
	512	rif = 0.0 ! flowste
	513	us = 0.0
	514	usws = 0.0
	515	vsws = 0.0
	516	ENDIF
	517
	518	ELSE
	519	e = 0.0 ! must be set, because used in
	520	rif = 0.0 ! flowste
	521	ts = 0.0
	522	us = 0.0
	523	usws = 0.0
	524	vsws = 0.0
	525	ENDIF
[102]	526	uswst = top_momentumflux_u
	527	vswst = top_momentumflux_v
[1]	528
	529	!
	530	!-- In every case qs = 0.0 (see also pt)
	531	!-- This could actually be computed more accurately in the 1D model.
	532	!-- Update when opportunity arises!
[75]	533	IF ( humidity .OR. passive_scalar ) qs = 0.0
[1]	534
	535	!
	536	!-- inside buildings set velocities back to zero
	537	IF ( topography /= 'flat' ) THEN
	538	DO i = nxl-1, nxr+1
	539	DO j = nys-1, nyn+1
	540	u(nzb:nzb_u_inner(j,i),j,i) = 0.0
	541	v(nzb:nzb_v_inner(j,i),j,i) = 0.0
	542	ENDDO
	543	ENDDO
[667]	544
[1]	545	!
	546	!-- WARNING: The extra boundary conditions set after running the
	547	!-- ------- 1D model impose an error on the divergence one layer
	548	!-- below the topography; need to correct later
	549	!-- ATTENTION: Provisional correction for Piacsek & Williams
	550	!-- --------- advection scheme: keep u and v zero one layer below
	551	!-- the topography.
	552	!
[667]	553	!-- Following was removed, because mirror boundary condition are
	554	!-- replaced by dirichlet boundary conditions
[1]	555	!
[667]	556	! IF ( ibc_uv_b == 0 ) THEN
	557	!!
	558	!!-- Satisfying the Dirichlet condition with an extra layer below
	559	!!-- the surface where the u and v component change their sign.
	560	! DO i = nxl-1, nxr+1
	561	! DO j = nys-1, nyn+1
	562	! IF ( nzb_u_inner(j,i) == 0 ) u(0,j,i) = -u(1,j,i)
	563	! IF ( nzb_v_inner(j,i) == 0 ) v(0,j,i) = -v(1,j,i)
	564	! ENDDO
	565	! ENDDO
	566	!
	567	! ELSE
	568	IF ( ibc_uv_b == 1 ) THEN
	569	!
[1]	570	!-- Neumann condition
	571	DO i = nxl-1, nxr+1
	572	DO j = nys-1, nyn+1
	573	IF ( nzb_u_inner(j,i) == 0 ) u(0,j,i) = u(1,j,i)
	574	IF ( nzb_v_inner(j,i) == 0 ) v(0,j,i) = v(1,j,i)
	575	ENDDO
	576	ENDDO
	577
	578	ENDIF
	579
	580	ENDIF
	581
	582	ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 ) &
	583	THEN
	584	!
	585	!-- Use constructed initial profiles (velocity constant with height,
	586	!-- temperature profile with constant gradient)
[667]	587	DO i = nxlg, nxrg
	588	DO j = nysg, nyng
[1]	589	pt(:,j,i) = pt_init
	590	u(:,j,i) = u_init
	591	v(:,j,i) = v_init
	592	ENDDO
	593	ENDDO
[75]	594
[1]	595	!
[292]	596	!-- Set initial horizontal velocities at the lowest computational grid
	597	!-- levels to zero in order to avoid too small time steps caused by the
	598	!-- diffusion limit in the initial phase of a run (at k=1, dz/2 occurs
	599	!-- in the limiting formula!). The original values are stored to be later
	600	!-- used for volume flow control.
[667]	601	DO i = nxlg, nxrg
	602	DO j = nysg, nyng
[1]	603	u(nzb:nzb_u_inner(j,i)+1,j,i) = 0.0
	604	v(nzb:nzb_v_inner(j,i)+1,j,i) = 0.0
	605	ENDDO
	606	ENDDO
	607
[75]	608	IF ( humidity .OR. passive_scalar ) THEN
[667]	609	DO i = nxlg, nxrg
	610	DO j = nysg, nyng
[1]	611	q(:,j,i) = q_init
	612	ENDDO
	613	ENDDO
	614	ENDIF
	615
[94]	616	IF ( ocean ) THEN
[667]	617	DO i = nxlg, nxrg
	618	DO j = nysg, nyng
[94]	619	sa(:,j,i) = sa_init
	620	ENDDO
	621	ENDDO
	622	ENDIF
[1]	623
	624	IF ( constant_diffusion ) THEN
	625	km = km_constant
	626	kh = km / prandtl_number
[108]	627	e = 0.0
	628	ELSEIF ( e_init > 0.0 ) THEN
	629	DO k = nzb+1, nzt
	630	km(k,:,:) = 0.1 * l_grid(k) * SQRT( e_init )
	631	ENDDO
	632	km(nzb,:,:) = km(nzb+1,:,:)
	633	km(nzt+1,:,:) = km(nzt,:,:)
	634	kh = km / prandtl_number
	635	e = e_init
[1]	636	ELSE
[108]	637	IF ( .NOT. ocean ) THEN
	638	kh = 0.01 ! there must exist an initial diffusion, because
	639	km = 0.01 ! otherwise no TKE would be produced by the
	640	! production terms, as long as not yet
	641	! e = (u/cm)*2 at k=nzb+1
	642	ELSE
	643	kh = 0.00001
	644	km = 0.00001
	645	ENDIF
	646	e = 0.0
[1]	647	ENDIF
[102]	648	rif = 0.0
	649	ts = 0.0
	650	us = 0.0
	651	usws = 0.0
	652	uswst = top_momentumflux_u
	653	vsws = 0.0
	654	vswst = top_momentumflux_v
[75]	655	IF ( humidity .OR. passive_scalar ) qs = 0.0
[1]	656
	657	!
	658	!-- Compute initial temperature field and other constants used in case
	659	!-- of a sloping surface
	660	IF ( sloping_surface ) CALL init_slope
	661
[46]	662	ELSEIF ( INDEX(initializing_actions, 'by_user') /= 0 ) &
	663	THEN
	664	!
	665	!-- Initialization will completely be done by the user
	666	CALL user_init_3d_model
	667
[1]	668	ENDIF
[667]	669	!
	670	!-- Bottom boundary
	671	IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2 ) THEN
	672	u(nzb,:,:) = 0.0
	673	v(nzb,:,:) = 0.0
	674	ENDIF
[1]	675
	676	!
[151]	677	!-- Apply channel flow boundary condition
[132]	678	IF ( TRIM( bc_uv_t ) == 'dirichlet_0' ) THEN
	679
	680	u(nzt+1,:,:) = 0.0
	681	v(nzt+1,:,:) = 0.0
	682
[151]	683	!-- For the Dirichlet condition to be correctly applied at the top, set
[132]	684	!-- ug and vg to zero there
	685	ug(nzt+1) = 0.0
	686	vg(nzt+1) = 0.0
	687
	688	ENDIF
	689
	690	!
[1]	691	!-- Calculate virtual potential temperature
[75]	692	IF ( humidity ) vpt = pt * ( 1.0 + 0.61 * q )
[1]	693
	694	!
	695	!-- Store initial profiles for output purposes etc.
	696	hom(:,1,5,:) = SPREAD( u(:,nys,nxl), 2, statistic_regions+1 )
	697	hom(:,1,6,:) = SPREAD( v(:,nys,nxl), 2, statistic_regions+1 )
[667]	698	IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2) THEN
	699	hom(nzb,1,5,:) = 0.0
	700	hom(nzb,1,6,:) = 0.0
[1]	701	ENDIF
	702	hom(:,1,7,:) = SPREAD( pt(:,nys,nxl), 2, statistic_regions+1 )
	703	hom(:,1,23,:) = SPREAD( km(:,nys,nxl), 2, statistic_regions+1 )
	704	hom(:,1,24,:) = SPREAD( kh(:,nys,nxl), 2, statistic_regions+1 )
	705
[97]	706	IF ( ocean ) THEN
	707	!
	708	!-- Store initial salinity profile
	709	hom(:,1,26,:) = SPREAD( sa(:,nys,nxl), 2, statistic_regions+1 )
	710	ENDIF
[1]	711
[75]	712	IF ( humidity ) THEN
[1]	713	!
	714	!-- Store initial profile of total water content, virtual potential
	715	!-- temperature
	716	hom(:,1,26,:) = SPREAD( q(:,nys,nxl), 2, statistic_regions+1 )
	717	hom(:,1,29,:) = SPREAD( vpt(:,nys,nxl), 2, statistic_regions+1 )
	718	IF ( cloud_physics .OR. cloud_droplets ) THEN
	719	!
	720	!-- Store initial profile of specific humidity and potential
	721	!-- temperature
	722	hom(:,1,27,:) = SPREAD( q(:,nys,nxl), 2, statistic_regions+1 )
	723	hom(:,1,28,:) = SPREAD( pt(:,nys,nxl), 2, statistic_regions+1 )
	724	ENDIF
	725	ENDIF
	726
	727	IF ( passive_scalar ) THEN
	728	!
	729	!-- Store initial scalar profile
	730	hom(:,1,26,:) = SPREAD( q(:,nys,nxl), 2, statistic_regions+1 )
	731	ENDIF
	732
	733	!
[19]	734	!-- Initialize fluxes at bottom surface
[1]	735	IF ( use_surface_fluxes ) THEN
	736
	737	IF ( constant_heatflux ) THEN
	738	!
	739	!-- Heat flux is prescribed
	740	IF ( random_heatflux ) THEN
	741	CALL disturb_heatflux
	742	ELSE
	743	shf = surface_heatflux
	744	!
	745	!-- Over topography surface_heatflux is replaced by wall_heatflux(0)
	746	IF ( TRIM( topography ) /= 'flat' ) THEN
[667]	747	DO i = nxlg, nxrg
	748	DO j = nysg, nyng
[1]	749	IF ( nzb_s_inner(j,i) /= 0 ) THEN
	750	shf(j,i) = wall_heatflux(0)
	751	ENDIF
	752	ENDDO
	753	ENDDO
	754	ENDIF
	755	ENDIF
	756	IF ( ASSOCIATED( shf_m ) ) shf_m = shf
	757	ENDIF
	758
	759	!
	760	!-- Determine the near-surface water flux
[75]	761	IF ( humidity .OR. passive_scalar ) THEN
[1]	762	IF ( constant_waterflux ) THEN
	763	qsws = surface_waterflux
[407]	764	!
	765	!-- Over topography surface_waterflux is replaced by
	766	!-- wall_humidityflux(0)
	767	IF ( TRIM( topography ) /= 'flat' ) THEN
	768	wall_qflux = wall_humidityflux
[667]	769	DO i = nxlg, nxrg
	770	DO j = nysg, nyng
[407]	771	IF ( nzb_s_inner(j,i) /= 0 ) THEN
	772	qsws(j,i) = wall_qflux(0)
	773	ENDIF
	774	ENDDO
	775	ENDDO
	776	ENDIF
[1]	777	IF ( ASSOCIATED( qsws_m ) ) qsws_m = qsws
	778	ENDIF
	779	ENDIF
	780
	781	ENDIF
	782
	783	!
[19]	784	!-- Initialize fluxes at top surface
[94]	785	!-- Currently, only the heatflux and salinity flux can be prescribed.
	786	!-- The latent flux is zero in this case!
[19]	787	IF ( use_top_fluxes ) THEN
	788
	789	IF ( constant_top_heatflux ) THEN
	790	!
	791	!-- Heat flux is prescribed
	792	tswst = top_heatflux
	793	IF ( ASSOCIATED( tswst_m ) ) tswst_m = tswst
	794
[75]	795	IF ( humidity .OR. passive_scalar ) THEN
[19]	796	qswst = 0.0
	797	IF ( ASSOCIATED( qswst_m ) ) qswst_m = qswst
	798	ENDIF
[94]	799
	800	IF ( ocean ) THEN
[95]	801	saswsb = bottom_salinityflux
[94]	802	saswst = top_salinityflux
	803	ENDIF
[102]	804	ENDIF
[19]	805
[102]	806	!
	807	!-- Initialization in case of a coupled model run
	808	IF ( coupling_mode == 'ocean_to_atmosphere' ) THEN
	809	tswst = 0.0
	810	IF ( ASSOCIATED( tswst_m ) ) tswst_m = tswst
	811	ENDIF
	812
[19]	813	ENDIF
	814
	815	!
[1]	816	!-- Initialize Prandtl layer quantities
	817	IF ( prandtl_layer ) THEN
	818
	819	z0 = roughness_length
	820
	821	IF ( .NOT. constant_heatflux ) THEN
	822	!
	823	!-- Surface temperature is prescribed. Here the heat flux cannot be
	824	!-- simply estimated, because therefore rif, u* and theta* would have
	825	!-- to be computed by iteration. This is why the heat flux is assumed
	826	!-- to be zero before the first time step. It approaches its correct
	827	!-- value in the course of the first few time steps.
	828	shf = 0.0
	829	IF ( ASSOCIATED( shf_m ) ) shf_m = 0.0
	830	ENDIF
	831
[75]	832	IF ( humidity .OR. passive_scalar ) THEN
[1]	833	IF ( .NOT. constant_waterflux ) THEN
	834	qsws = 0.0
	835	IF ( ASSOCIATED( qsws_m ) ) qsws_m = 0.0
	836	ENDIF
	837	ENDIF
	838
	839	ENDIF
	840
[152]	841
	842	!
[707]	843	!-- For the moment, vertical velocity is zero
	844	w = 0.0
[1]	845
	846	!
	847	!-- Initialize array sums (must be defined in first call of pres)
	848	sums = 0.0
	849
	850	!
[707]	851	!-- In case of iterative solvers, p must get an initial value
	852	IF ( psolver == 'multigrid' .OR. psolver == 'sor' ) p = 0.0
	853
	854	!
[72]	855	!-- Treating cloud physics, liquid water content and precipitation amount
	856	!-- are zero at beginning of the simulation
	857	IF ( cloud_physics ) THEN
	858	ql = 0.0
	859	IF ( precipitation ) precipitation_amount = 0.0
	860	ENDIF
[673]	861
	862	!
	863	!-- Initialize quantities for special advections schemes
	864	CALL init_advec
[1]	865
	866	!
	867	!-- Impose vortex with vertical axis on the initial velocity profile
	868	IF ( INDEX( initializing_actions, 'initialize_vortex' ) /= 0 ) THEN
	869	CALL init_rankine
	870	ENDIF
	871
	872	!
	873	!-- Impose temperature anomaly (advection test only)
	874	IF ( INDEX( initializing_actions, 'initialize_ptanom' ) /= 0 ) THEN
	875	CALL init_pt_anomaly
	876	ENDIF
	877
	878	!
	879	!-- If required, change the surface temperature at the start of the 3D run
	880	IF ( pt_surface_initial_change /= 0.0 ) THEN
	881	pt(nzb,:,:) = pt(nzb,:,:) + pt_surface_initial_change
	882	ENDIF
	883
	884	!
	885	!-- If required, change the surface humidity/scalar at the start of the 3D
	886	!-- run
[75]	887	IF ( ( humidity .OR. passive_scalar ) .AND. &
[1]	888	q_surface_initial_change /= 0.0 ) THEN
	889	q(nzb,:,:) = q(nzb,:,:) + q_surface_initial_change
	890	ENDIF
	891
	892	!
	893	!-- Initialize the random number generator (from numerical recipes)
	894	CALL random_function_ini
	895
	896	!
	897	!-- Initialize old and new time levels.
	898	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	899	e_m = e; pt_m = pt; u_m = u; v_m = v; w_m = w; kh_m = kh; km_m = km
	900	ELSE
	901	te_m = 0.0; tpt_m = 0.0; tu_m = 0.0; tv_m = 0.0; tw_m = 0.0
	902	ENDIF
	903	e_p = e; pt_p = pt; u_p = u; v_p = v; w_p = w
	904
[75]	905	IF ( humidity .OR. passive_scalar ) THEN
[1]	906	IF ( ASSOCIATED( q_m ) ) q_m = q
	907	IF ( timestep_scheme(1:5) == 'runge' ) tq_m = 0.0
	908	q_p = q
[75]	909	IF ( humidity .AND. ASSOCIATED( vpt_m ) ) vpt_m = vpt
[1]	910	ENDIF
	911
[94]	912	IF ( ocean ) THEN
	913	tsa_m = 0.0
	914	sa_p = sa
	915	ENDIF
[667]	916
[94]	917
[147]	918	ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data' .OR. &
[667]	919	TRIM( initializing_actions ) == 'cyclic_fill' ) &
[1]	920	THEN
	921	!
[328]	922	!-- When reading data for cyclic fill of 3D prerun data, first read
[147]	923	!-- some of the global variables from restart file
[328]	924	IF ( TRIM( initializing_actions ) == 'cyclic_fill' ) THEN
[559]	925
	926	CALL read_parts_of_var_list
[147]	927	CALL close_file( 13 )
[328]	928
[151]	929	!
[328]	930	!-- Initialization of the turbulence recycling method
	931	IF ( turbulent_inflow ) THEN
	932	!
	933	!-- Store temporally and horizontally averaged vertical profiles to be
	934	!-- used as mean inflow profiles
	935	ALLOCATE( mean_inflow_profiles(nzb:nzt+1,5) )
[151]	936
[328]	937	mean_inflow_profiles(:,1) = hom_sum(:,1,0) ! u
	938	mean_inflow_profiles(:,2) = hom_sum(:,2,0) ! v
	939	mean_inflow_profiles(:,4) = hom_sum(:,4,0) ! pt
	940	mean_inflow_profiles(:,5) = hom_sum(:,8,0) ! e
[151]	941
	942	!
[328]	943	!-- Use these mean profiles for the inflow (provided that Dirichlet
	944	!-- conditions are used)
	945	IF ( inflow_l ) THEN
[667]	946	DO j = nysg, nyng
[328]	947	DO k = nzb, nzt+1
[667]	948	u(k,j,nxlg:-1) = mean_inflow_profiles(k,1)
	949	v(k,j,nxlg:-1) = mean_inflow_profiles(k,2)
	950	w(k,j,nxlg:-1) = 0.0
	951	pt(k,j,nxlg:-1) = mean_inflow_profiles(k,4)
	952	e(k,j,nxlg:-1) = mean_inflow_profiles(k,5)
[328]	953	ENDDO
[151]	954	ENDDO
[328]	955	ENDIF
[151]	956
	957	!
[328]	958	!-- Calculate the damping factors to be used at the inflow. For a
	959	!-- turbulent inflow the turbulent fluctuations have to be limited
	960	!-- vertically because otherwise the turbulent inflow layer will grow
	961	!-- in time.
	962	IF ( inflow_damping_height == 9999999.9 ) THEN
[151]	963	!
[328]	964	!-- Default: use the inversion height calculated by the prerun; if
	965	!-- this is zero, inflow_damping_height must be explicitly
	966	!-- specified.
	967	IF ( hom_sum(nzb+6,pr_palm,0) /= 0.0 ) THEN
	968	inflow_damping_height = hom_sum(nzb+6,pr_palm,0)
	969	ELSE
	970	WRITE( message_string, * ) 'inflow_damping_height must be ',&
	971	'explicitly specified because&the inversion height ', &
	972	'calculated by the prerun is zero.'
	973	CALL message( 'init_3d_model', 'PA0318', 1, 2, 0, 6, 0 )
	974	ENDIF
	975
[292]	976	ENDIF
[151]	977
[328]	978	IF ( inflow_damping_width == 9999999.9 ) THEN
[151]	979	!
[328]	980	!-- Default for the transition range: one tenth of the undamped
	981	!-- layer
	982	inflow_damping_width = 0.1 * inflow_damping_height
[151]	983
[328]	984	ENDIF
[151]	985
[328]	986	ALLOCATE( inflow_damping_factor(nzb:nzt+1) )
[151]	987
[328]	988	DO k = nzb, nzt+1
[151]	989
[328]	990	IF ( zu(k) <= inflow_damping_height ) THEN
	991	inflow_damping_factor(k) = 1.0
	992	ELSEIF ( zu(k) <= inflow_damping_height + &
	993	inflow_damping_width ) THEN
	994	inflow_damping_factor(k) = 1.0 - &
[151]	995	( zu(k) - inflow_damping_height ) / &
	996	inflow_damping_width
[328]	997	ELSE
	998	inflow_damping_factor(k) = 0.0
	999	ENDIF
[151]	1000
[328]	1001	ENDDO
	1002	ENDIF
[151]	1003
[147]	1004	ENDIF
	1005
[152]	1006	!
[163]	1007	!-- Read binary data from restart file
[667]	1008
[559]	1009	CALL read_3d_binary
[163]	1010
	1011	!
[359]	1012	!-- Inside buildings set velocities and TKE back to zero
	1013	IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND. &
	1014	topography /= 'flat' ) THEN
	1015	!
	1016	!-- Inside buildings set velocities and TKE back to zero.
	1017	!-- Other scalars (pt, q, s, km, kh, p, sa, ...) are ignored at present,
	1018	!-- maybe revise later.
	1019	IF ( timestep_scheme(1:5) == 'runge' ) THEN
[667]	1020	DO i = nxlg, nxrg
	1021	DO j = nysg, nyng
[359]	1022	u (nzb:nzb_u_inner(j,i),j,i) = 0.0
	1023	v (nzb:nzb_v_inner(j,i),j,i) = 0.0
	1024	w (nzb:nzb_w_inner(j,i),j,i) = 0.0
	1025	e (nzb:nzb_w_inner(j,i),j,i) = 0.0
	1026	u_m(nzb:nzb_u_inner(j,i),j,i) = 0.0
	1027	v_m(nzb:nzb_v_inner(j,i),j,i) = 0.0
	1028	w_m(nzb:nzb_w_inner(j,i),j,i) = 0.0
	1029	e_m(nzb:nzb_w_inner(j,i),j,i) = 0.0
	1030	tu_m(nzb:nzb_u_inner(j,i),j,i) = 0.0
	1031	tv_m(nzb:nzb_v_inner(j,i),j,i) = 0.0
	1032	tw_m(nzb:nzb_w_inner(j,i),j,i) = 0.0
	1033	te_m(nzb:nzb_w_inner(j,i),j,i) = 0.0
	1034	tpt_m(nzb:nzb_w_inner(j,i),j,i) = 0.0
	1035	ENDDO
	1036	ENDDO
	1037	ELSE
[667]	1038	DO i = nxlg, nxrg
	1039	DO j = nysg, nyng
[359]	1040	u (nzb:nzb_u_inner(j,i),j,i) = 0.0
	1041	v (nzb:nzb_v_inner(j,i),j,i) = 0.0
	1042	w (nzb:nzb_w_inner(j,i),j,i) = 0.0
	1043	e (nzb:nzb_w_inner(j,i),j,i) = 0.0
	1044	u_m(nzb:nzb_u_inner(j,i),j,i) = 0.0
	1045	v_m(nzb:nzb_v_inner(j,i),j,i) = 0.0
	1046	w_m(nzb:nzb_w_inner(j,i),j,i) = 0.0
	1047	e_m(nzb:nzb_w_inner(j,i),j,i) = 0.0
	1048	u_p(nzb:nzb_u_inner(j,i),j,i) = 0.0
	1049	v_p(nzb:nzb_v_inner(j,i),j,i) = 0.0
	1050	w_p(nzb:nzb_w_inner(j,i),j,i) = 0.0
	1051	e_p(nzb:nzb_w_inner(j,i),j,i) = 0.0
	1052	ENDDO
	1053	ENDDO
	1054	ENDIF
	1055
	1056	ENDIF
	1057
	1058	!
[1]	1059	!-- Calculate initial temperature field and other constants used in case
	1060	!-- of a sloping surface
	1061	IF ( sloping_surface ) CALL init_slope
	1062
	1063	!
	1064	!-- Initialize new time levels (only done in order to set boundary values
	1065	!-- including ghost points)
	1066	e_p = e; pt_p = pt; u_p = u; v_p = v; w_p = w
[75]	1067	IF ( humidity .OR. passive_scalar ) q_p = q
[94]	1068	IF ( ocean ) sa_p = sa
[1]	1069
[181]	1070	!
	1071	!-- Allthough tendency arrays are set in prognostic_equations, they have
	1072	!-- have to be predefined here because they are used (but multiplied with 0)
	1073	!-- there before they are set.
	1074	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1075	te_m = 0.0; tpt_m = 0.0; tu_m = 0.0; tv_m = 0.0; tw_m = 0.0
	1076	IF ( humidity .OR. passive_scalar ) tq_m = 0.0
	1077	IF ( ocean ) tsa_m = 0.0
	1078	ENDIF
	1079
[1]	1080	ELSE
	1081	!
	1082	!-- Actually this part of the programm should not be reached
[254]	1083	message_string = 'unknown initializing problem'
	1084	CALL message( 'init_3d_model', 'PA0193', 1, 2, 0, 6, 0 )
[1]	1085	ENDIF
	1086
[151]	1087
	1088	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
[1]	1089	!
[151]	1090	!-- Initialize old timelevels needed for radiation boundary conditions
	1091	IF ( outflow_l ) THEN
	1092	u_m_l(:,:,:) = u(:,:,1:2)
	1093	v_m_l(:,:,:) = v(:,:,0:1)
	1094	w_m_l(:,:,:) = w(:,:,0:1)
	1095	ENDIF
	1096	IF ( outflow_r ) THEN
	1097	u_m_r(:,:,:) = u(:,:,nx-1:nx)
	1098	v_m_r(:,:,:) = v(:,:,nx-1:nx)
	1099	w_m_r(:,:,:) = w(:,:,nx-1:nx)
	1100	ENDIF
	1101	IF ( outflow_s ) THEN
	1102	u_m_s(:,:,:) = u(:,0:1,:)
	1103	v_m_s(:,:,:) = v(:,1:2,:)
	1104	w_m_s(:,:,:) = w(:,0:1,:)
	1105	ENDIF
	1106	IF ( outflow_n ) THEN
	1107	u_m_n(:,:,:) = u(:,ny-1:ny,:)
	1108	v_m_n(:,:,:) = v(:,ny-1:ny,:)
	1109	w_m_n(:,:,:) = w(:,ny-1:ny,:)
	1110	ENDIF
[667]	1111
[151]	1112	ENDIF
[680]	1113
[667]	1114	!
	1115	!-- Calculate the initial volume flow at the right and north boundary
	1116	IF ( conserve_volume_flow ) THEN
[151]	1117
[667]	1118	volume_flow_initial_l = 0.0
	1119	volume_flow_area_l = 0.0
	1120
	1121	IF ( TRIM( initializing_actions ) == 'cyclic_fill' ) THEN
	1122
	1123	IF ( nxr == nx ) THEN
	1124	DO j = nys, nyn
	1125	DO k = nzb_2d(j,nx) + 1, nzt
	1126	volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
	1127	hom_sum(k,1,0) * dzw(k)
	1128	volume_flow_area_l(1) = volume_flow_area_l(1) + dzw(k)
	1129	ENDDO
	1130	ENDDO
	1131	ENDIF
	1132
	1133	IF ( nyn == ny ) THEN
	1134	DO i = nxl, nxr
	1135	DO k = nzb_2d(ny,i) + 1, nzt
	1136	volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
	1137	hom_sum(k,2,0) * dzw(k)
	1138	volume_flow_area_l(2) = volume_flow_area_l(2) + dzw(k)
	1139	ENDDO
	1140	ENDDO
	1141	ENDIF
	1142
	1143	ELSEIF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
	1144
	1145	IF ( nxr == nx ) THEN
	1146	DO j = nys, nyn
	1147	DO k = nzb_2d(j,nx) + 1, nzt
	1148	volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
	1149	u(k,j,nx) * dzw(k)
	1150	volume_flow_area_l(1) = volume_flow_area_l(1) + dzw(k)
	1151	ENDDO
	1152	ENDDO
	1153	ENDIF
	1154
	1155	IF ( nyn == ny ) THEN
	1156	DO i = nxl, nxr
	1157	DO k = nzb_2d(ny,i) + 1, nzt
	1158	volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
	1159	v(k,ny,i) * dzw(k)
	1160	volume_flow_area_l(2) = volume_flow_area_l(2) + dzw(k)
	1161	ENDDO
	1162	ENDDO
	1163	ENDIF
	1164
	1165	ENDIF
	1166
	1167	#if defined( __parallel )
[680]	1168	CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&
	1169	2, MPI_REAL, MPI_SUM, comm2d, ierr )
	1170	CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1), &
	1171	2, MPI_REAL, MPI_SUM, comm2d, ierr )
[667]	1172
	1173	#else
[680]	1174	volume_flow_initial = volume_flow_initial_l
	1175	volume_flow_area = volume_flow_area_l
[667]	1176	#endif
	1177
[151]	1178	!
[667]	1179	!-- In case of 'bulk_velocity' mode, volume_flow_initial is overridden
	1180	!-- and calculated from u\|v_bulk instead.
[680]	1181	IF ( TRIM( conserve_volume_flow_mode ) == 'bulk_velocity' ) THEN
	1182	volume_flow_initial(1) = u_bulk * volume_flow_area(1)
	1183	volume_flow_initial(2) = v_bulk * volume_flow_area(2)
	1184	ENDIF
[667]	1185
[680]	1186	ENDIF
	1187
	1188
[667]	1189	!
[680]	1190	!-- Impose random perturbation on the horizontal velocity field and then
	1191	!-- remove the divergences from the velocity field at the initial stage
	1192	IF ( create_disturbances .AND. &
	1193	TRIM( initializing_actions ) /= 'read_restart_data' .AND. &
	1194	TRIM( initializing_actions ) /= 'cyclic_fill' ) THEN
	1195
	1196	CALL disturb_field( nzb_u_inner, tend, u )
	1197	CALL disturb_field( nzb_v_inner, tend, v )
	1198	n_sor = nsor_ini
	1199	CALL pres
	1200	n_sor = nsor
	1201	ENDIF
	1202
	1203	!
[138]	1204	!-- Initialization of the leaf area density
[680]	1205	IF ( plant_canopy ) THEN
[138]	1206
	1207	SELECT CASE ( TRIM( canopy_mode ) )
	1208
	1209	CASE( 'block' )
	1210
[667]	1211	DO i = nxlg, nxrg
	1212	DO j = nysg, nyng
[138]	1213	lad_s(:,j,i) = lad(:)
	1214	cdc(:,j,i) = drag_coefficient
[153]	1215	IF ( passive_scalar ) THEN
	1216	sls(:,j,i) = leaf_surface_concentration
	1217	sec(:,j,i) = scalar_exchange_coefficient
	1218	ENDIF
[138]	1219	ENDDO
	1220	ENDDO
	1221
	1222	CASE DEFAULT
	1223
	1224	!
	1225	!-- The DEFAULT case is reached either if the parameter
	1226	!-- canopy mode contains a wrong character string or if the
	1227	!-- user has coded a special case in the user interface.
	1228	!-- There, the subroutine user_init_plant_canopy checks
	1229	!-- which of these two conditions applies.
	1230	CALL user_init_plant_canopy
	1231
	1232	END SELECT
	1233
[667]	1234	CALL exchange_horiz( lad_s, nbgp )
	1235	CALL exchange_horiz( cdc, nbgp )
[138]	1236
[153]	1237	IF ( passive_scalar ) THEN
[667]	1238	CALL exchange_horiz( sls, nbgp )
	1239	CALL exchange_horiz( sec, nbgp )
[153]	1240	ENDIF
	1241
	1242	!
	1243	!-- Sharp boundaries of the plant canopy in horizontal directions
	1244	!-- In vertical direction the interpolation is retained, as the leaf
	1245	!-- area density is initialised by prescribing a vertical profile
	1246	!-- consisting of piecewise linear segments. The upper boundary
	1247	!-- of the plant canopy is now defined by lad_w(pch_index,:,:) = 0.0.
	1248
[138]	1249	DO i = nxl, nxr
	1250	DO j = nys, nyn
	1251	DO k = nzb, nzt+1
[153]	1252	IF ( lad_s(k,j,i) > 0.0 ) THEN
	1253	lad_u(k,j,i) = lad_s(k,j,i)
	1254	lad_u(k,j,i+1) = lad_s(k,j,i)
	1255	lad_v(k,j,i) = lad_s(k,j,i)
	1256	lad_v(k,j+1,i) = lad_s(k,j,i)
	1257	ENDIF
[138]	1258	ENDDO
	1259	DO k = nzb, nzt
	1260	lad_w(k,j,i) = 0.5 * ( lad_s(k+1,j,i) + lad_s(k,j,i) )
	1261	ENDDO
	1262	ENDDO
	1263	ENDDO
	1264
[153]	1265	lad_w(pch_index,:,:) = 0.0
	1266	lad_w(nzt+1,:,:) = lad_w(nzt,:,:)
[138]	1267
[667]	1268	CALL exchange_horiz( lad_u, nbgp )
	1269	CALL exchange_horiz( lad_v, nbgp )
	1270	CALL exchange_horiz( lad_w, nbgp )
[153]	1271
	1272	!
	1273	!-- Initialisation of the canopy heat source distribution
	1274	IF ( cthf /= 0.0 ) THEN
	1275	!
	1276	!-- Piecewise evaluation of the leaf area index by
	1277	!-- integration of the leaf area density
	1278	lai(:,:,:) = 0.0
[667]	1279	DO i = nxlg, nxrg
	1280	DO j = nysg, nyng
[153]	1281	DO k = pch_index-1, 0, -1
	1282	lai(k,j,i) = lai(k+1,j,i) + &
	1283	( 0.5 * ( lad_w(k+1,j,i) + &
	1284	lad_s(k+1,j,i) ) * &
	1285	( zw(k+1) - zu(k+1) ) ) + &
	1286	( 0.5 * ( lad_w(k,j,i) + &
	1287	lad_s(k+1,j,i) ) * &
	1288	( zu(k+1) - zw(k) ) )
	1289	ENDDO
	1290	ENDDO
	1291	ENDDO
	1292
	1293	!
	1294	!-- Evaluation of the upward kinematic vertical heat flux within the
	1295	!-- canopy
[667]	1296	DO i = nxlg, nxrg
	1297	DO j = nysg, nyng
[153]	1298	DO k = 0, pch_index
	1299	canopy_heat_flux(k,j,i) = cthf * &
	1300	exp( -0.6 * lai(k,j,i) )
	1301	ENDDO
	1302	ENDDO
	1303	ENDDO
	1304
	1305	!
	1306	!-- The near surface heat flux is derived from the heat flux
	1307	!-- distribution within the canopy
	1308	shf(:,:) = canopy_heat_flux(0,:,:)
	1309
	1310	IF ( ASSOCIATED( shf_m ) ) shf_m = shf
	1311
	1312	ENDIF
	1313
[138]	1314	ENDIF
	1315
	1316	!
[1]	1317	!-- If required, initialize dvrp-software
	1318	IF ( dt_dvrp /= 9999999.9 ) CALL init_dvrp
	1319
[96]	1320	IF ( ocean ) THEN
[1]	1321	!
[96]	1322	!-- Initialize quantities needed for the ocean model
	1323	CALL init_ocean
[388]	1324
[96]	1325	ELSE
	1326	!
	1327	!-- Initialize quantities for handling cloud physics
	1328	!-- This routine must be called before init_particles, because
	1329	!-- otherwise, array pt_d_t, needed in data_output_dvrp (called by
	1330	!-- init_particles) is not defined.
	1331	CALL init_cloud_physics
	1332	ENDIF
[1]	1333
	1334	!
	1335	!-- If required, initialize particles
[63]	1336	IF ( particle_advection ) CALL init_particles
[1]	1337
	1338	!
[673]	1339	!-- Initialize the ws-scheme.
	1340	IF ( ws_scheme_sca .OR. ws_scheme_mom ) CALL ws_init
[1]	1341
	1342	!
[673]	1343	!-- Setting weighting factors for right evaluation of perturbation pressure
	1344	!-- and turbulent quantities during the RK substeps.
	1345	IF ( TRIM(timestep_scheme) == 'runge-kutta-3' ) THEN ! RK3
	1346	weight_substep(1) = 0.166666666666666
	1347	weight_substep(2) = 0.3
	1348	weight_substep(3) = 0.533333333333333
	1349
	1350	weight_pres(1) = 0.333333333333333
	1351	weight_pres(2) = 0.416666666666666
	1352	weight_pres(3) = 0.25
	1353	ELSEIF ( TRIM(timestep_scheme) == 'runge-kutta-2' ) THEN ! RK2
	1354	weight_substep(1) = 0.5
	1355	weight_substep(2) = 0.5
	1356
	1357	weight_pres(1) = 0.5
	1358	weight_pres(2) = 0.5
	1359	ELSE ! Euler and Leapfrog
	1360	weight_substep(1) = 1.0
	1361	weight_pres(1) = 1.0
	1362	ENDIF
	1363
	1364	!
[1]	1365	!-- Initialize Rayleigh damping factors
	1366	rdf = 0.0
	1367	IF ( rayleigh_damping_factor /= 0.0 ) THEN
[108]	1368	IF ( .NOT. ocean ) THEN
	1369	DO k = nzb+1, nzt
	1370	IF ( zu(k) >= rayleigh_damping_height ) THEN
	1371	rdf(k) = rayleigh_damping_factor * &
[1]	1372	( SIN( pi * 0.5 * ( zu(k) - rayleigh_damping_height ) &
	1373	/ ( zu(nzt) - rayleigh_damping_height ) )&
	1374	)**2
[108]	1375	ENDIF
	1376	ENDDO
	1377	ELSE
	1378	DO k = nzt, nzb+1, -1
	1379	IF ( zu(k) <= rayleigh_damping_height ) THEN
	1380	rdf(k) = rayleigh_damping_factor * &
	1381	( SIN( pi * 0.5 * ( rayleigh_damping_height - zu(k) ) &
	1382	/ ( rayleigh_damping_height - zu(nzb+1)))&
	1383	)**2
	1384	ENDIF
	1385	ENDDO
	1386	ENDIF
[1]	1387	ENDIF
	1388
	1389	!
[240]	1390	!-- Initialize the starting level and the vertical smoothing factor used for
	1391	!-- the external pressure gradient
	1392	dp_smooth_factor = 1.0
	1393	IF ( dp_external ) THEN
	1394	!
	1395	!-- Set the starting level dp_level_ind_b only if it has not been set before
	1396	!-- (e.g. in init_grid).
	1397	IF ( dp_level_ind_b == 0 ) THEN
	1398	ind_array = MINLOC( ABS( dp_level_b - zu ) )
	1399	dp_level_ind_b = ind_array(1) - 1 + nzb
	1400	! MINLOC uses lower array bound 1
	1401	ENDIF
	1402	IF ( dp_smooth ) THEN
	1403	dp_smooth_factor(:dp_level_ind_b) = 0.0
	1404	DO k = dp_level_ind_b+1, nzt
	1405	dp_smooth_factor(k) = 0.5 * ( 1.0 + SIN( pi * &
	1406	( REAL( k - dp_level_ind_b ) / &
	1407	REAL( nzt - dp_level_ind_b ) - 0.5 ) ) )
	1408	ENDDO
	1409	ENDIF
	1410	ENDIF
	1411
	1412	!
[1]	1413	!-- Initialize diffusivities used within the outflow damping layer in case of
	1414	!-- non-cyclic lateral boundaries. A linear increase is assumed over the first
	1415	!-- half of the width of the damping layer
[707]	1416	IF ( bc_lr_dirrad ) THEN
[1]	1417
[667]	1418	DO i = nxlg, nxrg
[73]	1419	IF ( i >= nx - outflow_damping_width ) THEN
	1420	km_damp_x(i) = km_damp_max * MIN( 1.0, &
	1421	( i - ( nx - outflow_damping_width ) ) / &
	1422	REAL( outflow_damping_width/2 ) &
	1423	)
	1424	ELSE
	1425	km_damp_x(i) = 0.0
	1426	ENDIF
	1427	ENDDO
[1]	1428
[707]	1429	ELSEIF ( bc_lr_raddir ) THEN
[1]	1430
[667]	1431	DO i = nxlg, nxrg
[73]	1432	IF ( i <= outflow_damping_width ) THEN
	1433	km_damp_x(i) = km_damp_max * MIN( 1.0, &
	1434	( outflow_damping_width - i ) / &
	1435	REAL( outflow_damping_width/2 ) &
	1436	)
	1437	ELSE
	1438	km_damp_x(i) = 0.0
	1439	ENDIF
	1440	ENDDO
[1]	1441
[73]	1442	ENDIF
[1]	1443
[707]	1444	IF ( bc_ns_dirrad ) THEN
[1]	1445
[667]	1446	DO j = nysg, nyng
[73]	1447	IF ( j >= ny - outflow_damping_width ) THEN
	1448	km_damp_y(j) = km_damp_max * MIN( 1.0, &
	1449	( j - ( ny - outflow_damping_width ) ) / &
	1450	REAL( outflow_damping_width/2 ) &
	1451	)
	1452	ELSE
	1453	km_damp_y(j) = 0.0
[1]	1454	ENDIF
	1455	ENDDO
	1456
[707]	1457	ELSEIF ( bc_ns_raddir ) THEN
[1]	1458
[667]	1459	DO j = nysg, nyng
[73]	1460	IF ( j <= outflow_damping_width ) THEN
	1461	km_damp_y(j) = km_damp_max * MIN( 1.0, &
	1462	( outflow_damping_width - j ) / &
	1463	REAL( outflow_damping_width/2 ) &
	1464	)
	1465	ELSE
	1466	km_damp_y(j) = 0.0
[1]	1467	ENDIF
[73]	1468	ENDDO
[1]	1469
	1470	ENDIF
	1471
	1472	!
	1473	!-- Initialize local summation arrays for UP flow_statistics. This is necessary
	1474	!-- because they may not yet have been initialized when they are called from
	1475	!-- flow_statistics (or - depending on the chosen model run - are never
	1476	!-- initialized)
	1477	sums_divnew_l = 0.0
	1478	sums_divold_l = 0.0
	1479	sums_l_l = 0.0
	1480	sums_up_fraction_l = 0.0
	1481	sums_wsts_bc_l = 0.0
	1482
	1483	!
	1484	!-- Pre-set masks for regional statistics. Default is the total model domain.
	1485	rmask = 1.0
	1486
	1487	!
[51]	1488	!-- User-defined initializing actions. Check afterwards, if maximum number
	1489	!-- of allowed timeseries is not exceeded
[1]	1490	CALL user_init
	1491
[51]	1492	IF ( dots_num > dots_max ) THEN
[254]	1493	WRITE( message_string, * ) 'number of time series quantities exceeds', &
[274]	1494	' its maximum of dots_max = ', dots_max, &
[254]	1495	' &Please increase dots_max in modules.f90.'
	1496	CALL message( 'init_3d_model', 'PA0194', 1, 2, 0, 6, 0 )
[51]	1497	ENDIF
	1498
[1]	1499	!
	1500	!-- Input binary data file is not needed anymore. This line must be placed
	1501	!-- after call of user_init!
	1502	CALL close_file( 13 )
	1503
	1504	!
	1505	!-- Compute total sum of active mask grid points
	1506	!-- ngp_2dh: number of grid points of a horizontal cross section through the
	1507	!-- total domain
	1508	!-- ngp_3d: number of grid points of the total domain
[132]	1509	ngp_2dh_outer_l = 0
	1510	ngp_2dh_outer = 0
	1511	ngp_2dh_s_inner_l = 0
	1512	ngp_2dh_s_inner = 0
	1513	ngp_2dh_l = 0
	1514	ngp_2dh = 0
[485]	1515	ngp_3d_inner_l = 0.0
[132]	1516	ngp_3d_inner = 0
	1517	ngp_3d = 0
	1518	ngp_sums = ( nz + 2 ) * ( pr_palm + max_pr_user )
[1]	1519
	1520	DO sr = 0, statistic_regions
	1521	DO i = nxl, nxr
	1522	DO j = nys, nyn
	1523	IF ( rmask(j,i,sr) == 1.0 ) THEN
	1524	!
	1525	!-- All xy-grid points
	1526	ngp_2dh_l(sr) = ngp_2dh_l(sr) + 1
	1527	!
	1528	!-- xy-grid points above topography
	1529	DO k = nzb_s_outer(j,i), nz + 1
	1530	ngp_2dh_outer_l(k,sr) = ngp_2dh_outer_l(k,sr) + 1
	1531	ENDDO
[132]	1532	DO k = nzb_s_inner(j,i), nz + 1
	1533	ngp_2dh_s_inner_l(k,sr) = ngp_2dh_s_inner_l(k,sr) + 1
	1534	ENDDO
[1]	1535	!
	1536	!-- All grid points of the total domain above topography
	1537	ngp_3d_inner_l(sr) = ngp_3d_inner_l(sr) + &
	1538	( nz - nzb_s_inner(j,i) + 2 )
	1539	ENDIF
	1540	ENDDO
	1541	ENDDO
	1542	ENDDO
	1543
	1544	sr = statistic_regions + 1
	1545	#if defined( __parallel )
[622]	1546	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[485]	1547	CALL MPI_ALLREDUCE( ngp_2dh_l(0), ngp_2dh(0), sr, MPI_INTEGER, MPI_SUM, &
[1]	1548	comm2d, ierr )
[622]	1549	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[485]	1550	CALL MPI_ALLREDUCE( ngp_2dh_outer_l(0,0), ngp_2dh_outer(0,0), (nz+2)*sr, &
[1]	1551	MPI_INTEGER, MPI_SUM, comm2d, ierr )
[622]	1552	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[485]	1553	CALL MPI_ALLREDUCE( ngp_2dh_s_inner_l(0,0), ngp_2dh_s_inner(0,0), &
[132]	1554	(nz+2)*sr, MPI_INTEGER, MPI_SUM, comm2d, ierr )
[622]	1555	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[485]	1556	CALL MPI_ALLREDUCE( ngp_3d_inner_l(0), ngp_3d_inner_tmp(0), sr, MPI_REAL, &
[1]	1557	MPI_SUM, comm2d, ierr )
[485]	1558	ngp_3d_inner = INT( ngp_3d_inner_tmp, KIND = SELECTED_INT_KIND( 18 ) )
[1]	1559	#else
[132]	1560	ngp_2dh = ngp_2dh_l
	1561	ngp_2dh_outer = ngp_2dh_outer_l
	1562	ngp_2dh_s_inner = ngp_2dh_s_inner_l
[485]	1563	ngp_3d_inner = INT( ngp_3d_inner_l, KIND = SELECTED_INT_KIND( 18 ) )
[1]	1564	#endif
	1565
[560]	1566	ngp_3d = INT ( ngp_2dh, KIND = SELECTED_INT_KIND( 18 ) ) * &
	1567	INT ( (nz + 2 ), KIND = SELECTED_INT_KIND( 18 ) )
[1]	1568
	1569	!
	1570	!-- Set a lower limit of 1 in order to avoid zero divisions in flow_statistics,
	1571	!-- buoyancy, etc. A zero value will occur for cases where all grid points of
	1572	!-- the respective subdomain lie below the surface topography
[667]	1573	ngp_2dh_outer = MAX( 1, ngp_2dh_outer(:,:) )
[631]	1574	ngp_3d_inner = MAX( INT(1, KIND = SELECTED_INT_KIND( 18 )), &
	1575	ngp_3d_inner(:) )
[667]	1576	ngp_2dh_s_inner = MAX( 1, ngp_2dh_s_inner(:,:) )
[1]	1577
[485]	1578	DEALLOCATE( ngp_2dh_l, ngp_2dh_outer_l, ngp_3d_inner_l, ngp_3d_inner_tmp )
[1]	1579
	1580
	1581	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 @ 707

Download in other formats: