Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |