Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |