Home

Context Navigation

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

Last change on this file since 315 was 292, checked in by letzel, 16 years ago
Typographical errors (netcdf) If the inversion height calculated by the prerun is zero, inflow_damping_height must be explicitly specified (init_3d_model) Small bugfix concerning 3d 64bit netcdf output format (header) Adjustments for lcsgib and lcsgih (subjob)
Property svn:keywords set to `Id`
File size: 51.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |