Home

Context Navigation

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

Last change on this file since 1411 was 1411, checked in by suehring, 10 years ago
Initial horizontal velocity profiles were not set to zero at the first vertical
Property svn:keywords set to `Id`
File size: 67.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |