Home

Context Navigation

init_3d_model.f90 @ 1179

Last change on this file since 1179 was 1179, checked in by raasch, 11 years ago

New:
---
Initial profiles can be used as reference state in the buoyancy term. New parameter
reference_state introduced. Calculation and handling of reference state in buoyancy term revised.
binary version for restart files changed from 3.9 to 3.9a (no downward compatibility!),
initial profile for rho added to hom (id=77)

Errors:

small bugfix for background communication (time_integration)

Property svn:keywords set to Id

File size: 64.9 KB

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

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

Context Navigation

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

Download in other formats: