Home

Context Navigation

init_3d_model.f90 @ 4185

Last change on this file since 4185 was 4185, checked in by oliver.maas, 5 years ago

For initializing_actions = ' cyclic_fill': Overwrite u_init, v_init, pt_init, q_init and s_init with the (temporally) and horizontally averaged vertical profiles from the end of the prerun, because these profiles shall be used as the basic state for the rayleigh damping and the pt_damping.

Property svn:keywords set to Id

Property svn:mergeinfo set to False

/palm/branches/chemistry/SOURCE/init_3d_model.f90	2047-3190,3218-3297
/palm/branches/forwind/SOURCE/init_3d_model.f90	1564-1913
/palm/branches/mosaik_M2/init_3d_model.f90	2360-3471
/palm/branches/palm4u/SOURCE/init_3d_model.f90	2540-2692
/palm/branches/rans/SOURCE/init_3d_model.f90	2078-3128
/palm/branches/resler/SOURCE/init_3d_model.f90	2023-3605
/palm/branches/salsa/SOURCE/init_3d_model.f90	2503-3581

File size: 80.7 KB

Rev	Line
[1682]	1	!> @file init_3d_model.f90
[2000]	2	!------------------------------------------------------------------------------!
[2696]	3	! This file is part of the PALM model system.
[1036]	4	!
[2000]	5	! PALM is free software: you can redistribute it and/or modify it under the
	6	! terms of the GNU General Public License as published by the Free Software
	7	! Foundation, either version 3 of the License, or (at your option) any later
	8	! version.
[1036]	9	!
	10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
	11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
	12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
	13	!
	14	! You should have received a copy of the GNU General Public License along with
	15	! PALM. If not, see <http://www.gnu.org/licenses/>.
	16	!
[3648]	17	! Copyright 1997-2019 Leibniz Universitaet Hannover
[2000]	18	!------------------------------------------------------------------------------!
[1036]	19	!
[254]	20	! Current revisions:
[732]	21	! ------------------
[2233]	22	!
[3589]	23	!
[2233]	24	! Former revisions:
	25	! -----------------
	26	! $Id: init_3d_model.f90 4185 2019-08-23 13:49:38Z oliver.maas $
[4185]	27	! For initializing_actions = ' cyclic_fill':
	28	! Overwrite u_init, v_init, pt_init, q_init and s_init with the
	29	! (temporally) and horizontally averaged vertical profiles from the end
	30	! of the prerun, because these profiles shall be used as the basic state
	31	! for the rayleigh damping and the pt_damping.
	32	!
	33	! 4182 2019-08-22 15:20:23Z scharf
[4182]	34	! Corrected "Former revisions" section
	35	!
	36	! 4168 2019-08-16 13:50:17Z suehring
[4168]	37	! Replace function get_topography_top_index by topo_top_ind
	38	!
	39	! 4151 2019-08-09 08:24:30Z suehring
[4151]	40	! Add netcdf directive around input calls (fix for last commit)
	41	!
	42	! 4150 2019-08-08 20:00:47Z suehring
[4150]	43	! Input of additional surface variables independent on land- or urban-surface
	44	! model
	45	!
	46	! 4131 2019-08-02 11:06:18Z monakurppa
[4131]	47	! Allocate sums and sums_l to allow profile output for salsa variables.
	48	!
	49	! 4130 2019-08-01 13:04:13Z suehring
[4130]	50	! Effectively reduce 3D initialization to 1D initial profiles. This is because
	51	! 3D initialization produces structures in the w-component that are correlated
	52	! with the processor grid for some unknown reason
	53	!
	54	! 4090 2019-07-11 15:06:47Z Giersch
[4090]	55	! Unused variables removed
	56	!
	57	! 4088 2019-07-11 13:57:56Z Giersch
[4088]	58	! Pressure and density profile calculations revised using basic functions
	59	!
	60	! 4048 2019-06-21 21:00:21Z knoop
[4028]	61	! Further modularization of particle code components
	62	!
	63	! 4017 2019-06-06 12:16:46Z schwenkel
[3987]	64	! Convert most location messages to debug messages to reduce output in
	65	! job logfile to a minimum
	66	!
	67	!
[3939]	68	! unused variable removed
	69	!
	70	! 3937 2019-04-29 15:09:07Z suehring
[3937]	71	! Move initialization of synthetic turbulence generator behind initialization
	72	! of offline nesting. Remove call for stg_adjust, as this is now already done
	73	! in stg_init.
	74	!
	75	! 3900 2019-04-16 15:17:43Z suehring
[3900]	76	! Fix problem with LOD = 2 initialization
	77	!
	78	! 3885 2019-04-11 11:29:34Z kanani
[3885]	79	! Changes related to global restructuring of location messages and introduction
	80	! of additional debug messages
	81	!
	82	! 3849 2019-04-01 16:35:16Z knoop
[3747]	83	! Move initialization of rmask before initializing user_init_arrays
	84	!
	85	! 3711 2019-01-31 13:44:26Z knoop
[3711]	86	! Introduced module_interface_init_checks for post-init checks in modules
	87	!
	88	! 3700 2019-01-26 17:03:42Z knoop
[3685]	89	! Some interface calls moved to module_interface + cleanup
	90	!
	91	! 3648 2019-01-02 16:35:46Z suehring
[3648]	92	! Rename subroutines for surface-data output
[3569]	93	!
[4182]	94	! Revision 1.1 1998/03/09 16:22:22 raasch
	95	! Initial revision
	96	!
	97	!
[1]	98	! Description:
	99	! ------------
[1682]	100	!> Allocation of arrays and initialization of the 3D model via
	101	!> a) pre-run the 1D model
	102	!> or
	103	!> b) pre-set constant linear profiles
	104	!> or
	105	!> c) read values of a previous run
[1]	106	!------------------------------------------------------------------------------!
[1682]	107	SUBROUTINE init_3d_model
[1]	108
[3298]	109
[667]	110	USE advec_ws
[1320]	111
[1]	112	USE arrays_3d
[1849]	113
[3274]	114	USE basic_constants_and_equations_mod, &
[4090]	115	ONLY: c_p, g, l_v, pi, exner_function, exner_function_invers, &
[3274]	116	ideal_gas_law_rho, ideal_gas_law_rho_pt, barometric_formula
	117
	118	USE bulk_cloud_model_mod, &
[3685]	119	ONLY: bulk_cloud_model
[3274]	120
[3298]	121	USE chem_modules, &
[3685]	122	ONLY: max_pr_cs ! ToDo: this dependency needs to be removed cause it is ugly #new_dom
[3298]	123
[1]	124	USE control_parameters
[3298]	125
[1320]	126	USE grid_variables, &
[2037]	127	ONLY: dx, dy, ddx2_mg, ddy2_mg
[2817]	128
[1]	129	USE indices
[3469]	130
[1320]	131	USE kinds
[1496]	132
[2320]	133	USE lsf_nudging_mod, &
[3685]	134	ONLY: ls_forcing_surf
[1849]	135
[2338]	136	USE model_1d_mod, &
[3241]	137	ONLY: init_1d_model, l1d, u1d, v1d
[2338]	138
[3685]	139	USE module_interface, &
[3711]	140	ONLY: module_interface_init_arrays, &
	141	module_interface_init, &
	142	module_interface_init_checks
[3685]	143
[3159]	144	USE multi_agent_system_mod, &
	145	ONLY: agents_active, mas_init
	146
[1783]	147	USE netcdf_interface, &
[3700]	148	ONLY: dots_max
[2696]	149
[2906]	150	USE netcdf_data_input_mod, &
[4150]	151	ONLY: char_fill, &
	152	check_existence, &
	153	close_input_file, &
	154	get_attribute, &
	155	get_variable, &
	156	init_3d, &
	157	input_pids_static, &
	158	inquire_num_variables, &
	159	inquire_variable_names, &
	160	input_file_static, &
	161	netcdf_data_input_init_3d, &
	162	open_read_file, &
	163	real_2d
	164
[3347]	165	USE nesting_offl_mod, &
	166	ONLY: nesting_offl_init
[3294]	167
[1]	168	USE pegrid
[3298]	169
[3524]	170	#if defined( __parallel )
[2934]	171	USE pmc_interface, &
	172	ONLY: nested_run
[3524]	173	#endif
[2934]	174
[1320]	175	USE random_function_mod
[3685]	176
[1400]	177	USE random_generator_parallel, &
[2172]	178	ONLY: init_parallel_random_generator
[3685]	179
[2894]	180	USE read_restart_data_mod, &
[3685]	181	ONLY: rrd_read_parts_of_global, rrd_local
	182
[1320]	183	USE statistics, &
[1738]	184	ONLY: hom, hom_sum, mean_surface_level_height, pr_palm, rmask, &
[1833]	185	statistic_regions, sums, sums_divnew_l, sums_divold_l, sums_l, &
[2277]	186	sums_l_l, sums_wsts_bc_l, ts_value, &
[1833]	187	weight_pres, weight_substep
[2259]	188
	189	USE synthetic_turbulence_generator_mod, &
[3939]	190	ONLY: stg_init, use_syn_turb_gen
[3685]	191
[1691]	192	USE surface_layer_fluxes_mod, &
	193	ONLY: init_surface_layer_fluxes
[2232]	194
	195	USE surface_mod, &
[4150]	196	ONLY : init_single_surface_properties, &
	197	init_surface_arrays, &
	198	init_surfaces, &
	199	surf_def_h, &
	200	surf_def_v, &
	201	surf_lsm_h, &
[4168]	202	surf_usm_h
[3685]	203
[3849]	204	#if defined( _OPENACC )
	205	USE surface_mod, &
	206	ONLY : bc_h
	207	#endif
	208
[3648]	209	USE surface_data_output_mod, &
	210	ONLY: surface_data_output_init
[3685]	211
[2007]	212	USE transpose_indices
[1]	213
	214	IMPLICIT NONE
	215
[4150]	216	CHARACTER(LEN=50), DIMENSION(:), ALLOCATABLE :: vars_pids_static !< variable names in static input file
	217
	218	INTEGER(iwp) :: i !< grid index in x direction
	219	INTEGER(iwp) :: ind_array(1) !< dummy used to determine start index for external pressure forcing
	220	INTEGER(iwp) :: j !< grid index in y direction
	221	INTEGER(iwp) :: k !< grid index in z direction
	222	INTEGER(iwp) :: k_surf !< surface level index
	223	INTEGER(iwp) :: l !< running index over surface orientation
	224	INTEGER(iwp) :: m !< index of surface element in surface data type
	225	INTEGER(iwp) :: num_vars_pids_static !< number of variables in static input file
	226	INTEGER(iwp) :: nz_u_shift !< topography-top index on u-grid, used to vertically shift initial profiles
	227	INTEGER(iwp) :: nz_v_shift !< topography-top index on v-grid, used to vertically shift initial profiles
	228	INTEGER(iwp) :: nz_w_shift !< topography-top index on w-grid, used to vertically shift initial profiles
	229	INTEGER(iwp) :: nz_s_shift !< topography-top index on scalar-grid, used to vertically shift initial profiles
	230	INTEGER(iwp) :: nz_u_shift_l !< topography-top index on u-grid, used to vertically shift initial profiles
	231	INTEGER(iwp) :: nz_v_shift_l !< topography-top index on v-grid, used to vertically shift initial profiles
	232	INTEGER(iwp) :: nz_w_shift_l !< topography-top index on w-grid, used to vertically shift initial profiles
	233	INTEGER(iwp) :: nz_s_shift_l !< topography-top index on scalar-grid, used to vertically shift initial profiles
	234	INTEGER(iwp) :: nzt_l !< index of top PE boundary for multigrid level
	235	INTEGER(iwp) :: pids_static_id !< file id for static input file
	236	INTEGER(iwp) :: sr !< index of statistic region
[1]	237
[3547]	238	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: ngp_2dh_l !< toal number of horizontal grid points in statistical region on subdomain
[1]	239
[3547]	240	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: ngp_2dh_outer_l !< number of horizontal non-wall bounded grid points on subdomain
	241	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: ngp_2dh_s_inner_l !< number of horizontal non-topography grid points on subdomain
[1]	242
[4150]	243
	244
[3182]	245	REAL(wp), DIMENSION(:), ALLOCATABLE :: init_l !< dummy array used for averaging 3D data to obtain inital profiles
[2037]	246	REAL(wp), DIMENSION(:), ALLOCATABLE :: p_hydrostatic !< hydrostatic pressure
	247
	248	REAL(wp) :: dx_l !< grid spacing along x on different multigrid level
	249	REAL(wp) :: dy_l !< grid spacing along y on different multigrid level
	250
[3547]	251	REAL(wp), DIMENSION(1:3) :: volume_flow_area_l !< area of lateral and top model domain surface on local subdomain
	252	REAL(wp), DIMENSION(1:3) :: volume_flow_initial_l !< initial volume flow into model domain
[1]	253
[3547]	254	REAL(wp), DIMENSION(:), ALLOCATABLE :: mean_surface_level_height_l !< mean surface level height on subdomain
	255	REAL(wp), DIMENSION(:), ALLOCATABLE :: ngp_3d_inner_l !< total number of non-topography grid points on subdomain
	256	REAL(wp), DIMENSION(:), ALLOCATABLE :: ngp_3d_inner_tmp !< total number of non-topography grid points
[1]	257
[4150]	258	TYPE(real_2d) :: tmp_2d !< temporary variable to input additional surface-data from static file
	259
[3987]	260	CALL location_message( 'model initialization', 'start' )
	261
	262	IF ( debug_output ) CALL debug_message( 'allocating arrays', 'start' )
[1]	263	!
	264	!-- Allocate arrays
[1788]	265	ALLOCATE( mean_surface_level_height(0:statistic_regions), &
	266	mean_surface_level_height_l(0:statistic_regions), &
	267	ngp_2dh(0:statistic_regions), ngp_2dh_l(0:statistic_regions), &
	268	ngp_3d(0:statistic_regions), &
	269	ngp_3d_inner(0:statistic_regions), &
	270	ngp_3d_inner_l(0:statistic_regions), &
	271	ngp_3d_inner_tmp(0:statistic_regions), &
	272	sums_divnew_l(0:statistic_regions), &
[1]	273	sums_divold_l(0:statistic_regions) )
[1195]	274	ALLOCATE( dp_smooth_factor(nzb:nzt), rdf(nzb+1:nzt), rdf_sc(nzb+1:nzt) )
[1788]	275	ALLOCATE( ngp_2dh_outer(nzb:nzt+1,0:statistic_regions), &
	276	ngp_2dh_outer_l(nzb:nzt+1,0:statistic_regions), &
	277	ngp_2dh_s_inner(nzb:nzt+1,0:statistic_regions), &
	278	ngp_2dh_s_inner_l(nzb:nzt+1,0:statistic_regions), &
	279	rmask(nysg:nyng,nxlg:nxrg,0:statistic_regions), &
[4131]	280	sums(nzb:nzt+1,pr_palm+max_pr_user+max_pr_cs+max_pr_salsa), &
	281	sums_l(nzb:nzt+1,pr_palm+max_pr_user+max_pr_cs+max_pr_salsa,0:threads_per_task-1), &
[1788]	282	sums_l_l(nzb:nzt+1,0:statistic_regions,0:threads_per_task-1), &
[3700]	283	sums_wsts_bc_l(nzb:nzt+1,0:statistic_regions) )
	284	ALLOCATE( ts_value(dots_max,0:statistic_regions) )
[978]	285	ALLOCATE( ptdf_x(nxlg:nxrg), ptdf_y(nysg:nyng) )
[1]	286
[1788]	287	ALLOCATE( d(nzb+1:nzt,nys:nyn,nxl:nxr), &
	288	p(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
[1010]	289	tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	290
[2696]	291	ALLOCATE( pt_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
[1788]	292	pt_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	293	u_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	294	u_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	295	u_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	296	v_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	297	v_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	298	v_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	299	w_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	300	w_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
[667]	301	w_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1788]	302	IF ( .NOT. neutral ) THEN
[1032]	303	ALLOCATE( pt_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	304	ENDIF
[673]	305	!
[3747]	306	!-- Pre-set masks for regional statistics. Default is the total model domain.
	307	!-- Ghost points are excluded because counting values at the ghost boundaries
	308	!-- would bias the statistics
	309	rmask = 1.0_wp
	310	rmask(:,nxlg:nxl-1,:) = 0.0_wp; rmask(:,nxr+1:nxrg,:) = 0.0_wp
	311	rmask(nysg:nys-1,:,:) = 0.0_wp; rmask(nyn+1:nyng,:,:) = 0.0_wp
	312	!
[707]	313	!-- Following array is required for perturbation pressure within the iterative
	314	!-- pressure solvers. For the multistep schemes (Runge-Kutta), array p holds
	315	!-- the weighted average of the substeps and cannot be used in the Poisson
	316	!-- solver.
	317	IF ( psolver == 'sor' ) THEN
	318	ALLOCATE( p_loc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1575]	319	ELSEIF ( psolver(1:9) == 'multigrid' ) THEN
[707]	320	!
	321	!-- For performance reasons, multigrid is using one ghost layer only
	322	ALLOCATE( p_loc(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
[673]	323	ENDIF
[1]	324
[1111]	325	!
	326	!-- Array for storing constant coeffficients of the tridiagonal solver
	327	IF ( psolver == 'poisfft' ) THEN
[1212]	328	ALLOCATE( tri(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1,2) )
[1111]	329	ALLOCATE( tric(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1) )
	330	ENDIF
	331
[1960]	332	IF ( humidity ) THEN
[1]	333	!
[1960]	334	!-- 3D-humidity
[1788]	335	ALLOCATE( q_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	336	q_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
[3011]	337	q_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	338	vpt_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	339	ENDIF
[1960]	340
	341	IF ( passive_scalar ) THEN
[1]	342
[1960]	343	!
	344	!-- 3D scalar arrays
	345	ALLOCATE( s_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	346	s_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	347	s_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[3636]	348
[1960]	349	ENDIF
	350
[1]	351	!
[3302]	352	!-- Allocate and set 1d-profiles for Stokes drift velocity. It may be set to
	353	!-- non-zero values later in ocean_init
	354	ALLOCATE( u_stokes_zu(nzb:nzt+1), u_stokes_zw(nzb:nzt+1), &
	355	v_stokes_zu(nzb:nzt+1), v_stokes_zw(nzb:nzt+1) )
	356	u_stokes_zu(:) = 0.0_wp
	357	u_stokes_zw(:) = 0.0_wp
	358	v_stokes_zu(:) = 0.0_wp
	359	v_stokes_zw(:) = 0.0_wp
	360
	361	!
[2037]	362	!-- Allocation of anelastic and Boussinesq approximation specific arrays
	363	ALLOCATE( p_hydrostatic(nzb:nzt+1) )
	364	ALLOCATE( rho_air(nzb:nzt+1) )
	365	ALLOCATE( rho_air_zw(nzb:nzt+1) )
	366	ALLOCATE( drho_air(nzb:nzt+1) )
	367	ALLOCATE( drho_air_zw(nzb:nzt+1) )
	368	!
[4088]	369	!-- Density profile calculation for anelastic and Boussinesq approximation
	370	!-- In case of a Boussinesq approximation, a constant density is calculated
	371	!-- mainly for output purposes. This density do not need to be considered
	372	!-- in the model's system of equations.
	373	IF ( TRIM( approximation ) == 'anelastic' ) THEN
[2037]	374	DO k = nzb, nzt+1
[4088]	375	p_hydrostatic(k) = barometric_formula(zu(k), pt_surface * &
	376	exner_function(surface_pressure * 100.0_wp), &
	377	surface_pressure * 100.0_wp)
	378
	379	rho_air(k) = ideal_gas_law_rho_pt(p_hydrostatic(k), pt_init(k))
[2037]	380	ENDDO
[4088]	381
[2037]	382	DO k = nzb, nzt
	383	rho_air_zw(k) = 0.5_wp * ( rho_air(k) + rho_air(k+1) )
	384	ENDDO
[4088]	385
[2037]	386	rho_air_zw(nzt+1) = rho_air_zw(nzt) &
	387	+ 2.0_wp * ( rho_air(nzt+1) - rho_air_zw(nzt) )
[4088]	388
[2037]	389	ELSE
[2252]	390	DO k = nzb, nzt+1
[4088]	391	p_hydrostatic(k) = barometric_formula(zu(nzb), pt_surface * &
	392	exner_function(surface_pressure * 100.0_wp), &
	393	surface_pressure * 100.0_wp)
	394
	395	rho_air(k) = ideal_gas_law_rho_pt(p_hydrostatic(k), pt_init(nzb))
[2252]	396	ENDDO
[4088]	397
[2252]	398	DO k = nzb, nzt
	399	rho_air_zw(k) = 0.5_wp * ( rho_air(k) + rho_air(k+1) )
	400	ENDDO
[4088]	401
[2252]	402	rho_air_zw(nzt+1) = rho_air_zw(nzt) &
	403	+ 2.0_wp * ( rho_air(nzt+1) - rho_air_zw(nzt) )
[4088]	404
[2037]	405	ENDIF
[2696]	406	!
[2037]	407	!-- compute the inverse density array in order to avoid expencive divisions
	408	drho_air = 1.0_wp / rho_air
	409	drho_air_zw = 1.0_wp / rho_air_zw
	410
	411	!
	412	!-- Allocation of flux conversion arrays
	413	ALLOCATE( heatflux_input_conversion(nzb:nzt+1) )
	414	ALLOCATE( waterflux_input_conversion(nzb:nzt+1) )
	415	ALLOCATE( momentumflux_input_conversion(nzb:nzt+1) )
	416	ALLOCATE( heatflux_output_conversion(nzb:nzt+1) )
	417	ALLOCATE( waterflux_output_conversion(nzb:nzt+1) )
	418	ALLOCATE( momentumflux_output_conversion(nzb:nzt+1) )
	419
	420	!
	421	!-- calculate flux conversion factors according to approximation and in-/output mode
	422	DO k = nzb, nzt+1
	423
	424	IF ( TRIM( flux_input_mode ) == 'kinematic' ) THEN
	425	heatflux_input_conversion(k) = rho_air_zw(k)
	426	waterflux_input_conversion(k) = rho_air_zw(k)
	427	momentumflux_input_conversion(k) = rho_air_zw(k)
	428	ELSEIF ( TRIM( flux_input_mode ) == 'dynamic' ) THEN
[3274]	429	heatflux_input_conversion(k) = 1.0_wp / c_p
[2037]	430	waterflux_input_conversion(k) = 1.0_wp / l_v
	431	momentumflux_input_conversion(k) = 1.0_wp
	432	ENDIF
	433
	434	IF ( TRIM( flux_output_mode ) == 'kinematic' ) THEN
	435	heatflux_output_conversion(k) = drho_air_zw(k)
	436	waterflux_output_conversion(k) = drho_air_zw(k)
	437	momentumflux_output_conversion(k) = drho_air_zw(k)
	438	ELSEIF ( TRIM( flux_output_mode ) == 'dynamic' ) THEN
[3274]	439	heatflux_output_conversion(k) = c_p
[2037]	440	waterflux_output_conversion(k) = l_v
	441	momentumflux_output_conversion(k) = 1.0_wp
	442	ENDIF
	443
	444	IF ( .NOT. humidity ) THEN
	445	waterflux_input_conversion(k) = 1.0_wp
	446	waterflux_output_conversion(k) = 1.0_wp
	447	ENDIF
	448
	449	ENDDO
	450
	451	!
	452	!-- In case of multigrid method, compute grid lengths and grid factors for the
	453	!-- grid levels with respective density on each grid
	454	IF ( psolver(1:9) == 'multigrid' ) THEN
	455
	456	ALLOCATE( ddx2_mg(maximum_grid_level) )
	457	ALLOCATE( ddy2_mg(maximum_grid_level) )
	458	ALLOCATE( dzu_mg(nzb+1:nzt+1,maximum_grid_level) )
	459	ALLOCATE( dzw_mg(nzb+1:nzt+1,maximum_grid_level) )
	460	ALLOCATE( f1_mg(nzb+1:nzt,maximum_grid_level) )
	461	ALLOCATE( f2_mg(nzb+1:nzt,maximum_grid_level) )
	462	ALLOCATE( f3_mg(nzb+1:nzt,maximum_grid_level) )
	463	ALLOCATE( rho_air_mg(nzb:nzt+1,maximum_grid_level) )
	464	ALLOCATE( rho_air_zw_mg(nzb:nzt+1,maximum_grid_level) )
	465
	466	dzu_mg(:,maximum_grid_level) = dzu
	467	rho_air_mg(:,maximum_grid_level) = rho_air
	468	!
	469	!-- Next line to ensure an equally spaced grid.
	470	dzu_mg(1,maximum_grid_level) = dzu(2)
	471	rho_air_mg(nzb,maximum_grid_level) = rho_air(nzb) + &
	472	(rho_air(nzb) - rho_air(nzb+1))
	473
	474	dzw_mg(:,maximum_grid_level) = dzw
	475	rho_air_zw_mg(:,maximum_grid_level) = rho_air_zw
	476	nzt_l = nzt
	477	DO l = maximum_grid_level-1, 1, -1
	478	dzu_mg(nzb+1,l) = 2.0_wp * dzu_mg(nzb+1,l+1)
	479	dzw_mg(nzb+1,l) = 2.0_wp * dzw_mg(nzb+1,l+1)
	480	rho_air_mg(nzb,l) = rho_air_mg(nzb,l+1) + (rho_air_mg(nzb,l+1) - rho_air_mg(nzb+1,l+1))
	481	rho_air_zw_mg(nzb,l) = rho_air_zw_mg(nzb,l+1) + (rho_air_zw_mg(nzb,l+1) - rho_air_zw_mg(nzb+1,l+1))
	482	rho_air_mg(nzb+1,l) = rho_air_mg(nzb+1,l+1)
	483	rho_air_zw_mg(nzb+1,l) = rho_air_zw_mg(nzb+1,l+1)
	484	nzt_l = nzt_l / 2
	485	DO k = 2, nzt_l+1
	486	dzu_mg(k,l) = dzu_mg(2k-2,l+1) + dzu_mg(2k-1,l+1)
	487	dzw_mg(k,l) = dzw_mg(2k-2,l+1) + dzw_mg(2k-1,l+1)
	488	rho_air_mg(k,l) = rho_air_mg(2*k-1,l+1)
	489	rho_air_zw_mg(k,l) = rho_air_zw_mg(2*k-1,l+1)
	490	ENDDO
	491	ENDDO
	492
	493	nzt_l = nzt
	494	dx_l = dx
	495	dy_l = dy
	496	DO l = maximum_grid_level, 1, -1
	497	ddx2_mg(l) = 1.0_wp / dx_l**2
	498	ddy2_mg(l) = 1.0_wp / dy_l**2
	499	DO k = nzb+1, nzt_l
	500	f2_mg(k,l) = rho_air_zw_mg(k,l) / ( dzu_mg(k+1,l) * dzw_mg(k,l) )
	501	f3_mg(k,l) = rho_air_zw_mg(k-1,l) / ( dzu_mg(k,l) * dzw_mg(k,l) )
	502	f1_mg(k,l) = 2.0_wp * ( ddx2_mg(l) + ddy2_mg(l) ) &
	503	* rho_air_mg(k,l) + f2_mg(k,l) + f3_mg(k,l)
	504	ENDDO
	505	nzt_l = nzt_l / 2
	506	dx_l = dx_l * 2.0_wp
	507	dy_l = dy_l * 2.0_wp
	508	ENDDO
	509
	510	ENDIF
	511
	512	!
[1299]	513	!-- 1D-array for large scale subsidence velocity
[1361]	514	IF ( .NOT. ALLOCATED( w_subs ) ) THEN
	515	ALLOCATE ( w_subs(nzb:nzt+1) )
	516	w_subs = 0.0_wp
	517	ENDIF
[1299]	518
	519	!
[106]	520	!-- Arrays to store velocity data from t-dt and the phase speeds which
	521	!-- are needed for radiation boundary conditions
[3182]	522	IF ( bc_radiation_l ) THEN
[1788]	523	ALLOCATE( u_m_l(nzb:nzt+1,nysg:nyng,1:2), &
	524	v_m_l(nzb:nzt+1,nysg:nyng,0:1), &
[667]	525	w_m_l(nzb:nzt+1,nysg:nyng,0:1) )
[73]	526	ENDIF
[3182]	527	IF ( bc_radiation_r ) THEN
[1788]	528	ALLOCATE( u_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx), &
	529	v_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx), &
[667]	530	w_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx) )
[73]	531	ENDIF
[3182]	532	IF ( bc_radiation_l .OR. bc_radiation_r ) THEN
[1788]	533	ALLOCATE( c_u(nzb:nzt+1,nysg:nyng), c_v(nzb:nzt+1,nysg:nyng), &
[667]	534	c_w(nzb:nzt+1,nysg:nyng) )
[106]	535	ENDIF
[3182]	536	IF ( bc_radiation_s ) THEN
[1788]	537	ALLOCATE( u_m_s(nzb:nzt+1,0:1,nxlg:nxrg), &
	538	v_m_s(nzb:nzt+1,1:2,nxlg:nxrg), &
[667]	539	w_m_s(nzb:nzt+1,0:1,nxlg:nxrg) )
[73]	540	ENDIF
[3182]	541	IF ( bc_radiation_n ) THEN
[1788]	542	ALLOCATE( u_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg), &
	543	v_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg), &
[667]	544	w_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg) )
[73]	545	ENDIF
[3182]	546	IF ( bc_radiation_s .OR. bc_radiation_n ) THEN
[1788]	547	ALLOCATE( c_u(nzb:nzt+1,nxlg:nxrg), c_v(nzb:nzt+1,nxlg:nxrg), &
[667]	548	c_w(nzb:nzt+1,nxlg:nxrg) )
[106]	549	ENDIF
[3182]	550	IF ( bc_radiation_l .OR. bc_radiation_r .OR. bc_radiation_s .OR. &
	551	bc_radiation_n ) THEN
[978]	552	ALLOCATE( c_u_m_l(nzb:nzt+1), c_v_m_l(nzb:nzt+1), c_w_m_l(nzb:nzt+1) )
	553	ALLOCATE( c_u_m(nzb:nzt+1), c_v_m(nzb:nzt+1), c_w_m(nzb:nzt+1) )
	554	ENDIF
[73]	555
	556	!
[1]	557	!-- Initial assignment of the pointers
[1032]	558	IF ( .NOT. neutral ) THEN
	559	pt => pt_1; pt_p => pt_2; tpt_m => pt_3
	560	ELSE
	561	pt => pt_1; pt_p => pt_1; tpt_m => pt_3
	562	ENDIF
[1001]	563	u => u_1; u_p => u_2; tu_m => u_3
	564	v => v_1; v_p => v_2; tv_m => v_3
	565	w => w_1; w_p => w_2; tw_m => w_3
[1]	566
[1960]	567	IF ( humidity ) THEN
[1001]	568	q => q_1; q_p => q_2; tq_m => q_3
[3274]	569	vpt => vpt_1
[1001]	570	ENDIF
[1960]	571
	572	IF ( passive_scalar ) THEN
	573	s => s_1; s_p => s_2; ts_m => s_3
	574	ENDIF
[1]	575
	576	!
[2696]	577	!-- Initialize surface arrays
[2232]	578	CALL init_surface_arrays
	579	!
[3294]	580	!-- Allocate arrays for other modules
[3685]	581	CALL module_interface_init_arrays
[1551]	582
[1914]	583
[2320]	584	!
[709]	585	!-- Allocate arrays containing the RK coefficient for calculation of
	586	!-- perturbation pressure and turbulent fluxes. At this point values are
	587	!-- set for pressure calculation during initialization (where no timestep
	588	!-- is done). Further below the values needed within the timestep scheme
	589	!-- will be set.
[1788]	590	ALLOCATE( weight_substep(1:intermediate_timestep_count_max), &
[1878]	591	weight_pres(1:intermediate_timestep_count_max) )
[1340]	592	weight_substep = 1.0_wp
	593	weight_pres = 1.0_wp
[1918]	594	intermediate_timestep_count = 0 ! needed when simulated_time = 0.0
[673]	595
[3987]	596	IF ( debug_output ) CALL debug_message( 'allocating arrays', 'end' )
[1918]	597
[673]	598	!
[3014]	599	!-- Initialize time series
	600	ts_value = 0.0_wp
	601
	602	!
[1918]	603	!-- Initialize local summation arrays for routine flow_statistics.
	604	!-- This is necessary because they may not yet have been initialized when they
	605	!-- are called from flow_statistics (or - depending on the chosen model run -
	606	!-- are never initialized)
	607	sums_divnew_l = 0.0_wp
	608	sums_divold_l = 0.0_wp
	609	sums_l_l = 0.0_wp
	610	sums_wsts_bc_l = 0.0_wp
[3182]	611
[1918]	612	!
[1]	613	!-- Initialize model variables
[1788]	614	IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. &
[328]	615	TRIM( initializing_actions ) /= 'cyclic_fill' ) THEN
[1]	616	!
[2696]	617	!-- Initialization with provided input data derived from larger-scale model
	618	IF ( INDEX( initializing_actions, 'inifor' ) /= 0 ) THEN
[3987]	619	IF ( debug_output ) CALL debug_message( 'initializing with INIFOR', 'start' )
[2696]	620	!
[3051]	621	!-- Read initial 1D profiles or 3D data from NetCDF file, depending
	622	!-- on the provided level-of-detail.
[2696]	623	!-- At the moment, only u, v, w, pt and q are provided.
	624	CALL netcdf_data_input_init_3d
	625	!
[3182]	626	!-- Please note, Inifor provides data from nzb+1 to nzt.
	627	!-- Bottom and top boundary conditions for Inifor profiles are already
	628	!-- set (just after reading), so that this is not necessary here.
	629	!-- Depending on the provided level-of-detail, initial Inifor data is
	630	!-- either stored on data type (lod=1), or directly on 3D arrays (lod=2).
	631	!-- In order to obtain also initial profiles in case of lod=2 (which
	632	!-- is required for e.g. damping), average over 3D data.
	633	IF( init_3d%lod_u == 1 ) THEN
	634	u_init = init_3d%u_init
	635	ELSEIF( init_3d%lod_u == 2 ) THEN
	636	ALLOCATE( init_l(nzb:nzt+1) )
	637	DO k = nzb, nzt+1
	638	init_l(k) = SUM( u(k,nys:nyn,nxl:nxr) )
	639	ENDDO
	640	init_l = init_l / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
[1384]	641
[3182]	642	#if defined( __parallel )
	643	CALL MPI_ALLREDUCE( init_l, u_init, nzt+1-nzb+1, &
	644	MPI_REAL, MPI_SUM, comm2d, ierr )
	645	#else
	646	u_init = init_l
	647	#endif
	648	DEALLOCATE( init_l )
[3051]	649
[2696]	650	ENDIF
[3182]	651
	652	IF( init_3d%lod_v == 1 ) THEN
	653	v_init = init_3d%v_init
	654	ELSEIF( init_3d%lod_v == 2 ) THEN
	655	ALLOCATE( init_l(nzb:nzt+1) )
	656	DO k = nzb, nzt+1
	657	init_l(k) = SUM( v(k,nys:nyn,nxl:nxr) )
	658	ENDDO
	659	init_l = init_l / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
[2696]	660
[3182]	661	#if defined( __parallel )
	662	CALL MPI_ALLREDUCE( init_l, v_init, nzt+1-nzb+1, &
	663	MPI_REAL, MPI_SUM, comm2d, ierr )
	664	#else
	665	v_init = init_l
	666	#endif
	667	DEALLOCATE( init_l )
	668	ENDIF
	669	IF( .NOT. neutral ) THEN
	670	IF( init_3d%lod_pt == 1 ) THEN
	671	pt_init = init_3d%pt_init
	672	ELSEIF( init_3d%lod_pt == 2 ) THEN
	673	ALLOCATE( init_l(nzb:nzt+1) )
	674	DO k = nzb, nzt+1
	675	init_l(k) = SUM( pt(k,nys:nyn,nxl:nxr) )
	676	ENDDO
	677	init_l = init_l / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
	678
	679	#if defined( __parallel )
	680	CALL MPI_ALLREDUCE( init_l, pt_init, nzt+1-nzb+1, &
	681	MPI_REAL, MPI_SUM, comm2d, ierr )
	682	#else
	683	pt_init = init_l
	684	#endif
	685	DEALLOCATE( init_l )
	686	ENDIF
	687	ENDIF
	688
	689
	690	IF( humidity ) THEN
	691	IF( init_3d%lod_q == 1 ) THEN
	692	q_init = init_3d%q_init
	693	ELSEIF( init_3d%lod_q == 2 ) THEN
	694	ALLOCATE( init_l(nzb:nzt+1) )
	695	DO k = nzb, nzt+1
	696	init_l(k) = SUM( q(k,nys:nyn,nxl:nxr) )
	697	ENDDO
	698	init_l = init_l / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
	699
	700	#if defined( __parallel )
	701	CALL MPI_ALLREDUCE( init_l, q_init, nzt+1-nzb+1, &
	702	MPI_REAL, MPI_SUM, comm2d, ierr )
	703	#else
	704	q_init = init_l
	705	#endif
	706	DEALLOCATE( init_l )
	707	ENDIF
	708	ENDIF
	709
[2696]	710	!
[4130]	711	!-- Write initial profiles onto 3D arrays.
	712	!-- Work-around, 3D initialization of u,v,w creates artificial
	713	!-- structures wich correlate with the processor grid. The reason
	714	!-- for this is still unknown. To work-around this, 3D initialization
	715	!-- will be effectively reduce to a 1D initialization where no such
	716	!-- artificial structures appear.
[2696]	717	DO i = nxlg, nxrg
	718	DO j = nysg, nyng
[4130]	719	IF( init_3d%lod_u == 1 .OR. init_3d%lod_u == 2 ) &
	720	u(:,j,i) = u_init(:)
	721	IF( init_3d%lod_v == 1 .OR. init_3d%lod_u == 2 ) &
	722	v(:,j,i) = v_init(:)
	723	IF( .NOT. neutral .AND. &
	724	( init_3d%lod_pt == 1 .OR. init_3d%lod_pt == 2 ) ) &
[3051]	725	pt(:,j,i) = pt_init(:)
[4130]	726	IF( humidity .AND. &
	727	( init_3d%lod_q == 1 .OR. init_3d%lod_q == 2 ) ) &
	728	q(:,j,i) = q_init(:)
[2696]	729	ENDDO
	730	ENDDO
	731	!
[3182]	732	!-- Set geostrophic wind components.
[2938]	733	IF ( init_3d%from_file_ug ) THEN
	734	ug(:) = init_3d%ug_init(:)
	735	ENDIF
	736	IF ( init_3d%from_file_vg ) THEN
	737	vg(:) = init_3d%vg_init(:)
	738	ENDIF
[3404]	739	!
	740	!-- Set bottom and top boundary condition for geostrophic wind
[2938]	741	ug(nzt+1) = ug(nzt)
	742	vg(nzt+1) = vg(nzt)
[3404]	743	ug(nzb) = ug(nzb+1)
	744	vg(nzb) = vg(nzb+1)
[2696]	745	!
	746	!-- Set inital w to 0
	747	w = 0.0_wp
	748
	749	IF ( passive_scalar ) THEN
	750	DO i = nxlg, nxrg
	751	DO j = nysg, nyng
	752	s(:,j,i) = s_init
	753	ENDDO
	754	ENDDO
	755	ENDIF
	756
	757	!
	758	!-- Set velocity components at non-atmospheric / oceanic grid points to
	759	!-- zero.
	760	u = MERGE( u, 0.0_wp, BTEST( wall_flags_0, 1 ) )
	761	v = MERGE( v, 0.0_wp, BTEST( wall_flags_0, 2 ) )
	762	w = MERGE( w, 0.0_wp, BTEST( wall_flags_0, 3 ) )
[2700]	763	!
	764	!-- Initialize surface variables, e.g. friction velocity, momentum
	765	!-- fluxes, etc.
	766	CALL init_surfaces
[2696]	767
[3987]	768	IF ( debug_output ) CALL debug_message( 'initializing with INIFOR', 'end' )
[2696]	769	!
	770	!-- Initialization via computed 1D-model profiles
	771	ELSEIF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 ) THEN
	772
[3987]	773	IF ( debug_output ) CALL debug_message( 'initializing with 1D model profiles', 'start' )
[1]	774	!
	775	!-- Use solutions of the 1D model as initial profiles,
	776	!-- start 1D model
	777	CALL init_1d_model
	778	!
	779	!-- Transfer initial profiles to the arrays of the 3D model
[667]	780	DO i = nxlg, nxrg
	781	DO j = nysg, nyng
[1]	782	pt(:,j,i) = pt_init
	783	u(:,j,i) = u1d
	784	v(:,j,i) = v1d
	785	ENDDO
	786	ENDDO
	787
[1960]	788	IF ( humidity ) THEN
[667]	789	DO i = nxlg, nxrg
	790	DO j = nysg, nyng
[1]	791	q(:,j,i) = q_init
	792	ENDDO
	793	ENDDO
	794	ENDIF
[2292]	795
[1960]	796	IF ( passive_scalar ) THEN
	797	DO i = nxlg, nxrg
	798	DO j = nysg, nyng
	799	s(:,j,i) = s_init
	800	ENDDO
	801	ENDDO
	802	ENDIF
[1]	803	!
	804	!-- Store initial profiles for output purposes etc.
[2696]	805	IF ( .NOT. constant_diffusion ) THEN
[1]	806	hom(:,1,25,:) = SPREAD( l1d, 2, statistic_regions+1 )
	807	ENDIF
	808	!
[2696]	809	!-- Set velocities back to zero
[2758]	810	u = MERGE( u, 0.0_wp, BTEST( wall_flags_0, 1 ) )
	811	v = MERGE( v, 0.0_wp, BTEST( wall_flags_0, 2 ) )
[1]	812	!
[2696]	813	!-- WARNING: The extra boundary conditions set after running the
	814	!-- ------- 1D model impose an error on the divergence one layer
	815	!-- below the topography; need to correct later
	816	!-- ATTENTION: Provisional correction for Piacsek & Williams
	817	!-- --------- advection scheme: keep u and v zero one layer below
	818	!-- the topography.
	819	IF ( ibc_uv_b == 1 ) THEN
[667]	820	!
[2696]	821	!-- Neumann condition
	822	DO i = nxl-1, nxr+1
	823	DO j = nys-1, nyn+1
	824	u(nzb,j,i) = u(nzb+1,j,i)
	825	v(nzb,j,i) = v(nzb+1,j,i)
[1]	826	ENDDO
[2696]	827	ENDDO
[1]	828
	829	ENDIF
[2618]	830	!
	831	!-- Initialize surface variables, e.g. friction velocity, momentum
	832	!-- fluxes, etc.
	833	CALL init_surfaces
[1]	834
[3987]	835	IF ( debug_output ) CALL debug_message( 'initializing with 1D model profiles', 'end' )
[1384]	836
[1788]	837	ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 ) &
[1]	838	THEN
[1241]	839
[3987]	840	IF ( debug_output ) CALL debug_message( 'initializing with constant profiles', 'start' )
[2259]	841
	842	!
[1]	843	!-- Use constructed initial profiles (velocity constant with height,
	844	!-- temperature profile with constant gradient)
[667]	845	DO i = nxlg, nxrg
	846	DO j = nysg, nyng
[1]	847	pt(:,j,i) = pt_init
	848	u(:,j,i) = u_init
	849	v(:,j,i) = v_init
	850	ENDDO
	851	ENDDO
	852	!
[2758]	853	!-- Mask topography
	854	u = MERGE( u, 0.0_wp, BTEST( wall_flags_0, 1 ) )
	855	v = MERGE( v, 0.0_wp, BTEST( wall_flags_0, 2 ) )
	856	!
[292]	857	!-- Set initial horizontal velocities at the lowest computational grid
	858	!-- levels to zero in order to avoid too small time steps caused by the
	859	!-- diffusion limit in the initial phase of a run (at k=1, dz/2 occurs
[2758]	860	!-- in the limiting formula!).
	861	!-- Please note, in case land- or urban-surface model is used and a
	862	!-- spinup is applied, masking the lowest computational level is not
	863	!-- possible as MOST as well as energy-balance parametrizations will not
	864	!-- work with zero wind velocity.
	865	IF ( ibc_uv_b /= 1 .AND. .NOT. spinup ) THEN
[1815]	866	DO i = nxlg, nxrg
	867	DO j = nysg, nyng
[2232]	868	DO k = nzb, nzt
	869	u(k,j,i) = MERGE( u(k,j,i), 0.0_wp, &
	870	BTEST( wall_flags_0(k,j,i), 20 ) )
	871	v(k,j,i) = MERGE( v(k,j,i), 0.0_wp, &
	872	BTEST( wall_flags_0(k,j,i), 21 ) )
	873	ENDDO
[1815]	874	ENDDO
	875	ENDDO
	876	ENDIF
[1]	877
[1960]	878	IF ( humidity ) THEN
[667]	879	DO i = nxlg, nxrg
	880	DO j = nysg, nyng
[1]	881	q(:,j,i) = q_init
	882	ENDDO
	883	ENDDO
	884	ENDIF
[1960]	885
	886	IF ( passive_scalar ) THEN
	887	DO i = nxlg, nxrg
	888	DO j = nysg, nyng
	889	s(:,j,i) = s_init
	890	ENDDO
	891	ENDDO
	892	ENDIF
[1]	893
[1920]	894	!
[1]	895	!-- Compute initial temperature field and other constants used in case
	896	!-- of a sloping surface
	897	IF ( sloping_surface ) CALL init_slope
[2618]	898	!
	899	!-- Initialize surface variables, e.g. friction velocity, momentum
	900	!-- fluxes, etc.
	901	CALL init_surfaces
[3579]	902
[3987]	903	IF ( debug_output ) CALL debug_message( 'initializing with constant profiles', 'end' )
[1384]	904
[1788]	905	ELSEIF ( INDEX(initializing_actions, 'by_user') /= 0 ) &
[46]	906	THEN
[1384]	907
[3987]	908	IF ( debug_output ) CALL debug_message( 'initializing by user', 'start' )
[46]	909	!
[2618]	910	!-- Pre-initialize surface variables, i.e. setting start- and end-indices
	911	!-- at each (j,i)-location. Please note, this does not supersede
	912	!-- user-defined initialization of surface quantities.
	913	CALL init_surfaces
	914	!
[46]	915	!-- Initialization will completely be done by the user
	916	CALL user_init_3d_model
	917
[3987]	918	IF ( debug_output ) CALL debug_message( 'initializing by user', 'end' )
[1384]	919
[1]	920	ENDIF
[1384]	921
[3987]	922	IF ( debug_output ) CALL debug_message( 'initializing statistics, boundary conditions, etc.', 'start' )
[1384]	923
[667]	924	!
	925	!-- Bottom boundary
	926	IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2 ) THEN
[1340]	927	u(nzb,:,:) = 0.0_wp
	928	v(nzb,:,:) = 0.0_wp
[667]	929	ENDIF
[1]	930
	931	!
[151]	932	!-- Apply channel flow boundary condition
[132]	933	IF ( TRIM( bc_uv_t ) == 'dirichlet_0' ) THEN
[1340]	934	u(nzt+1,:,:) = 0.0_wp
	935	v(nzt+1,:,:) = 0.0_wp
[132]	936	ENDIF
	937
	938	!
[1]	939	!-- Calculate virtual potential temperature
[1960]	940	IF ( humidity ) vpt = pt * ( 1.0_wp + 0.61_wp * q )
[1]	941
	942	!
[2696]	943	!-- Store initial profiles for output purposes etc.. Please note, in case of
	944	!-- initialization of u, v, w, pt, and q via output data derived from larger
	945	!-- scale models, data will not be horizontally homogeneous. Actually, a mean
	946	!-- profile should be calculated before.
[1]	947	hom(:,1,5,:) = SPREAD( u(:,nys,nxl), 2, statistic_regions+1 )
	948	hom(:,1,6,:) = SPREAD( v(:,nys,nxl), 2, statistic_regions+1 )
[667]	949	IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2) THEN
[1340]	950	hom(nzb,1,5,:) = 0.0_wp
	951	hom(nzb,1,6,:) = 0.0_wp
[1]	952	ENDIF
	953	hom(:,1,7,:) = SPREAD( pt(:,nys,nxl), 2, statistic_regions+1 )
	954
[75]	955	IF ( humidity ) THEN
[1]	956	!
	957	!-- Store initial profile of total water content, virtual potential
	958	!-- temperature
	959	hom(:,1,26,:) = SPREAD( q(:,nys,nxl), 2, statistic_regions+1 )
	960	hom(:,1,29,:) = SPREAD( vpt(:,nys,nxl), 2, statistic_regions+1 )
[2696]	961	!
[3040]	962	!-- Store initial profile of mixing ratio and potential
[2696]	963	!-- temperature
[3274]	964	IF ( bulk_cloud_model .OR. cloud_droplets ) THEN
[1]	965	hom(:,1,27,:) = SPREAD( q(:,nys,nxl), 2, statistic_regions+1 )
	966	hom(:,1,28,:) = SPREAD( pt(:,nys,nxl), 2, statistic_regions+1 )
	967	ENDIF
	968	ENDIF
	969
[2696]	970	!
	971	!-- Store initial scalar profile
[1]	972	IF ( passive_scalar ) THEN
[2513]	973	hom(:,1,121,:) = SPREAD( s(:,nys,nxl), 2, statistic_regions+1 )
[1]	974	ENDIF
	975
	976	!
[1400]	977	!-- Initialize the random number generators (from numerical recipes)
	978	CALL random_function_ini
[1429]	979
[1400]	980	IF ( random_generator == 'random-parallel' ) THEN
[3241]	981	CALL init_parallel_random_generator( nx, nys, nyn, nxl, nxr )
[1400]	982	ENDIF
	983	!
[1179]	984	!-- Set the reference state to be used in the buoyancy terms (for ocean runs
	985	!-- the reference state will be set (overwritten) in init_ocean)
	986	IF ( use_single_reference_value ) THEN
[1788]	987	IF ( .NOT. humidity ) THEN
[1179]	988	ref_state(:) = pt_reference
	989	ELSE
	990	ref_state(:) = vpt_reference
	991	ENDIF
	992	ELSE
[1788]	993	IF ( .NOT. humidity ) THEN
[1179]	994	ref_state(:) = pt_init(:)
	995	ELSE
	996	ref_state(:) = vpt(:,nys,nxl)
	997	ENDIF
	998	ENDIF
[152]	999
	1000	!
[707]	1001	!-- For the moment, vertical velocity is zero
[1340]	1002	w = 0.0_wp
[1]	1003
	1004	!
	1005	!-- Initialize array sums (must be defined in first call of pres)
[1340]	1006	sums = 0.0_wp
[1]	1007
	1008	!
[707]	1009	!-- In case of iterative solvers, p must get an initial value
[1575]	1010	IF ( psolver(1:9) == 'multigrid' .OR. psolver == 'sor' ) p = 0.0_wp
[707]	1011	!
[1]	1012	!-- Impose vortex with vertical axis on the initial velocity profile
	1013	IF ( INDEX( initializing_actions, 'initialize_vortex' ) /= 0 ) THEN
	1014	CALL init_rankine
	1015	ENDIF
	1016
	1017	!
[3035]	1018	!-- Impose temperature anomaly (advection test only) or warm air bubble
	1019	!-- close to surface
	1020	IF ( INDEX( initializing_actions, 'initialize_ptanom' ) /= 0 .OR. &
	1021	INDEX( initializing_actions, 'initialize_bubble' ) /= 0 ) THEN
[1]	1022	CALL init_pt_anomaly
	1023	ENDIF
[3035]	1024
[1]	1025	!
	1026	!-- If required, change the surface temperature at the start of the 3D run
[1340]	1027	IF ( pt_surface_initial_change /= 0.0_wp ) THEN
[1]	1028	pt(nzb,:,:) = pt(nzb,:,:) + pt_surface_initial_change
	1029	ENDIF
	1030
	1031	!
	1032	!-- If required, change the surface humidity/scalar at the start of the 3D
	1033	!-- run
[1960]	1034	IF ( humidity .AND. q_surface_initial_change /= 0.0_wp ) &
[1]	1035	q(nzb,:,:) = q(nzb,:,:) + q_surface_initial_change
[1960]	1036
	1037	IF ( passive_scalar .AND. s_surface_initial_change /= 0.0_wp ) &
	1038	s(nzb,:,:) = s(nzb,:,:) + s_surface_initial_change
	1039
[1]	1040
	1041	!
	1042	!-- Initialize old and new time levels.
[2696]	1043	tpt_m = 0.0_wp; tu_m = 0.0_wp; tv_m = 0.0_wp; tw_m = 0.0_wp
	1044	pt_p = pt; u_p = u; v_p = v; w_p = w
[1]	1045
[1960]	1046	IF ( humidity ) THEN
[1340]	1047	tq_m = 0.0_wp
[1]	1048	q_p = q
	1049	ENDIF
[1960]	1050
	1051	IF ( passive_scalar ) THEN
	1052	ts_m = 0.0_wp
	1053	s_p = s
	1054	ENDIF
[1]	1055
[3987]	1056	IF ( debug_output ) CALL debug_message( 'initializing statistics, boundary conditions, etc.', 'end' )
[94]	1057
[1788]	1058	ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data' .OR. &
[2232]	1059	TRIM( initializing_actions ) == 'cyclic_fill' ) &
[1]	1060	THEN
[1384]	1061
[3987]	1062	IF ( debug_output ) CALL debug_message( 'initializing in case of restart / cyclic_fill', 'start' )
[1]	1063	!
[3609]	1064	!-- Initialize surface elements and its attributes, e.g. heat- and
	1065	!-- momentumfluxes, roughness, scaling parameters. As number of surface
	1066	!-- elements might be different between runs, e.g. in case of cyclic fill,
	1067	!-- and not all surface elements are read, surface elements need to be
	1068	!-- initialized before.
	1069	!-- Please note, in case of cyclic fill, surfaces should be initialized
	1070	!-- after restart data is read, else, individual settings of surface
	1071	!-- parameters will be overwritten from data of precursor run, hence,
	1072	!-- init_surfaces is called a second time after reading the restart data.
	1073	CALL init_surfaces
	1074	!
[767]	1075	!-- When reading data for cyclic fill of 3D prerun data files, read
	1076	!-- some of the global variables from the restart file which are required
	1077	!-- for initializing the inflow
[328]	1078	IF ( TRIM( initializing_actions ) == 'cyclic_fill' ) THEN
[559]	1079
[759]	1080	DO i = 0, io_blocks-1
	1081	IF ( i == io_group ) THEN
[2894]	1082	CALL rrd_read_parts_of_global
[759]	1083	ENDIF
	1084	#if defined( __parallel )
	1085	CALL MPI_BARRIER( comm2d, ierr )
	1086	#endif
	1087	ENDDO
[328]	1088
[767]	1089	ENDIF
	1090
[151]	1091	!
[2894]	1092	!-- Read processor specific binary data from restart file
[767]	1093	DO i = 0, io_blocks-1
	1094	IF ( i == io_group ) THEN
[2894]	1095	CALL rrd_local
[767]	1096	ENDIF
	1097	#if defined( __parallel )
	1098	CALL MPI_BARRIER( comm2d, ierr )
	1099	#endif
	1100	ENDDO
[3608]	1101	!
[3609]	1102	!-- In case of cyclic fill, call init_surfaces a second time, so that
	1103	!-- surface properties such as heat fluxes are initialized as prescribed.
	1104	IF ( TRIM( initializing_actions ) == 'cyclic_fill' ) &
	1105	CALL init_surfaces
[767]	1106
[328]	1107	!
[2550]	1108	!-- In case of complex terrain and cyclic fill method as initialization,
	1109	!-- shift initial data in the vertical direction for each point in the
	1110	!-- x-y-plane depending on local surface height
	1111	IF ( complex_terrain .AND. &
	1112	TRIM( initializing_actions ) == 'cyclic_fill' ) THEN
	1113	DO i = nxlg, nxrg
	1114	DO j = nysg, nyng
[4168]	1115	nz_u_shift = topo_top_ind(j,i,1)
	1116	nz_v_shift = topo_top_ind(j,i,2)
	1117	nz_w_shift = topo_top_ind(j,i,3)
	1118	nz_s_shift = topo_top_ind(j,i,0)
[2550]	1119
	1120	u(nz_u_shift:nzt+1,j,i) = u(0:nzt+1-nz_u_shift,j,i)
	1121
	1122	v(nz_v_shift:nzt+1,j,i) = v(0:nzt+1-nz_v_shift,j,i)
	1123
	1124	w(nz_w_shift:nzt+1,j,i) = w(0:nzt+1-nz_w_shift,j,i)
	1125
	1126	p(nz_s_shift:nzt+1,j,i) = p(0:nzt+1-nz_s_shift,j,i)
	1127	pt(nz_s_shift:nzt+1,j,i) = pt(0:nzt+1-nz_s_shift,j,i)
	1128	ENDDO
	1129	ENDDO
	1130	ENDIF
	1131	!
[767]	1132	!-- Initialization of the turbulence recycling method
[1788]	1133	IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND. &
[767]	1134	turbulent_inflow ) THEN
	1135	!
	1136	!-- First store the profiles to be used at the inflow.
	1137	!-- These profiles are the (temporally) and horizontally averaged vertical
	1138	!-- profiles from the prerun. Alternatively, prescribed profiles
[4185]	1139	!-- for u,v-components can be used.
	1140	!-- Overwrite u_init, v_init, pt_init, q_init and s_init with the
	1141	!-- (temporally) and horizontally averaged vertical profiles from the end
	1142	!-- of the prerun, because these profiles shall be used as the basic state
	1143	!-- for the rayleigh damping and the pt_damping.
[3288]	1144	ALLOCATE( mean_inflow_profiles(nzb:nzt+1,1:num_mean_inflow_profiles) )
[151]	1145
[767]	1146	IF ( use_prescribed_profile_data ) THEN
	1147	mean_inflow_profiles(:,1) = u_init ! u
	1148	mean_inflow_profiles(:,2) = v_init ! v
	1149	ELSE
[328]	1150	mean_inflow_profiles(:,1) = hom_sum(:,1,0) ! u
[4185]	1151	u_init(:) = hom_sum(:,1,0)
[328]	1152	mean_inflow_profiles(:,2) = hom_sum(:,2,0) ! v
[4185]	1153	v_init(:) = hom_sum(:,2,0)
[767]	1154	ENDIF
	1155	mean_inflow_profiles(:,4) = hom_sum(:,4,0) ! pt
[4185]	1156	pt_init(:) = hom_sum(:,4,0)
[1960]	1157	IF ( humidity ) &
	1158	mean_inflow_profiles(:,6) = hom_sum(:,41,0) ! q
[4185]	1159	q_init(:) = hom_sum(:,41,0)
[1960]	1160	IF ( passive_scalar ) &
	1161	mean_inflow_profiles(:,7) = hom_sum(:,115,0) ! s
[4185]	1162	s_init(:) = hom_sum(:,115,0)
[2550]	1163	!
	1164	!-- In case of complex terrain, determine vertical displacement at inflow
	1165	!-- boundary and adjust mean inflow profiles
	1166	IF ( complex_terrain ) THEN
	1167	IF ( nxlg <= 0 .AND. nxrg >= 0 .AND. nysg <= 0 .AND. nyng >= 0 ) THEN
[4168]	1168	nz_u_shift_l = topo_top_ind(j,i,1)
	1169	nz_v_shift_l = topo_top_ind(j,i,2)
	1170	nz_w_shift_l = topo_top_ind(j,i,3)
	1171	nz_s_shift_l = topo_top_ind(j,i,0)
[2550]	1172	ELSE
	1173	nz_u_shift_l = 0
	1174	nz_v_shift_l = 0
	1175	nz_w_shift_l = 0
	1176	nz_s_shift_l = 0
	1177	ENDIF
[151]	1178
[2550]	1179	#if defined( __parallel )
	1180	CALL MPI_ALLREDUCE(nz_u_shift_l, nz_u_shift, 1, MPI_INTEGER, &
	1181	MPI_MAX, comm2d, ierr)
	1182	CALL MPI_ALLREDUCE(nz_v_shift_l, nz_v_shift, 1, MPI_INTEGER, &
	1183	MPI_MAX, comm2d, ierr)
	1184	CALL MPI_ALLREDUCE(nz_w_shift_l, nz_w_shift, 1, MPI_INTEGER, &
	1185	MPI_MAX, comm2d, ierr)
	1186	CALL MPI_ALLREDUCE(nz_s_shift_l, nz_s_shift, 1, MPI_INTEGER, &
	1187	MPI_MAX, comm2d, ierr)
	1188	#else
	1189	nz_u_shift = nz_u_shift_l
	1190	nz_v_shift = nz_v_shift_l
	1191	nz_w_shift = nz_w_shift_l
	1192	nz_s_shift = nz_s_shift_l
	1193	#endif
	1194
	1195	mean_inflow_profiles(:,1) = 0.0_wp
	1196	mean_inflow_profiles(nz_u_shift:nzt+1,1) = hom_sum(0:nzt+1-nz_u_shift,1,0) ! u
	1197
	1198	mean_inflow_profiles(:,2) = 0.0_wp
	1199	mean_inflow_profiles(nz_v_shift:nzt+1,2) = hom_sum(0:nzt+1-nz_v_shift,2,0) ! v
	1200
	1201	mean_inflow_profiles(nz_s_shift:nzt+1,4) = hom_sum(0:nzt+1-nz_s_shift,4,0) ! pt
	1202
	1203	ENDIF
	1204
[151]	1205	!
[767]	1206	!-- If necessary, adjust the horizontal flow field to the prescribed
	1207	!-- profiles
	1208	IF ( use_prescribed_profile_data ) THEN
	1209	DO i = nxlg, nxrg
[667]	1210	DO j = nysg, nyng
[328]	1211	DO k = nzb, nzt+1
[767]	1212	u(k,j,i) = u(k,j,i) - hom_sum(k,1,0) + u_init(k)
	1213	v(k,j,i) = v(k,j,i) - hom_sum(k,2,0) + v_init(k)
[328]	1214	ENDDO
[151]	1215	ENDDO
[767]	1216	ENDDO
	1217	ENDIF
[151]	1218
	1219	!
[767]	1220	!-- Use these mean profiles at the inflow (provided that Dirichlet
	1221	!-- conditions are used)
[3182]	1222	IF ( bc_dirichlet_l ) THEN
[767]	1223	DO j = nysg, nyng
	1224	DO k = nzb, nzt+1
	1225	u(k,j,nxlg:-1) = mean_inflow_profiles(k,1)
	1226	v(k,j,nxlg:-1) = mean_inflow_profiles(k,2)
[1340]	1227	w(k,j,nxlg:-1) = 0.0_wp
[767]	1228	pt(k,j,nxlg:-1) = mean_inflow_profiles(k,4)
[1960]	1229	IF ( humidity ) &
[1615]	1230	q(k,j,nxlg:-1) = mean_inflow_profiles(k,6)
[1960]	1231	IF ( passive_scalar ) &
	1232	s(k,j,nxlg:-1) = mean_inflow_profiles(k,7)
[767]	1233	ENDDO
	1234	ENDDO
	1235	ENDIF
	1236
[151]	1237	!
[767]	1238	!-- Calculate the damping factors to be used at the inflow. For a
	1239	!-- turbulent inflow the turbulent fluctuations have to be limited
	1240	!-- vertically because otherwise the turbulent inflow layer will grow
	1241	!-- in time.
[1340]	1242	IF ( inflow_damping_height == 9999999.9_wp ) THEN
[767]	1243	!
	1244	!-- Default: use the inversion height calculated by the prerun; if
	1245	!-- this is zero, inflow_damping_height must be explicitly
	1246	!-- specified.
[1340]	1247	IF ( hom_sum(nzb+6,pr_palm,0) /= 0.0_wp ) THEN
[767]	1248	inflow_damping_height = hom_sum(nzb+6,pr_palm,0)
	1249	ELSE
[1788]	1250	WRITE( message_string, * ) 'inflow_damping_height must be ', &
	1251	'explicitly specified because&the inversion height ', &
[767]	1252	'calculated by the prerun is zero.'
	1253	CALL message( 'init_3d_model', 'PA0318', 1, 2, 0, 6, 0 )
[292]	1254	ENDIF
[151]	1255
[767]	1256	ENDIF
	1257
[1340]	1258	IF ( inflow_damping_width == 9999999.9_wp ) THEN
[151]	1259	!
[767]	1260	!-- Default for the transition range: one tenth of the undamped
	1261	!-- layer
[1340]	1262	inflow_damping_width = 0.1_wp * inflow_damping_height
[151]	1263
[767]	1264	ENDIF
[151]	1265
[767]	1266	ALLOCATE( inflow_damping_factor(nzb:nzt+1) )
[151]	1267
[767]	1268	DO k = nzb, nzt+1
[151]	1269
[767]	1270	IF ( zu(k) <= inflow_damping_height ) THEN
[1340]	1271	inflow_damping_factor(k) = 1.0_wp
[996]	1272	ELSEIF ( zu(k) <= ( inflow_damping_height + inflow_damping_width ) ) THEN
[1340]	1273	inflow_damping_factor(k) = 1.0_wp - &
[996]	1274	( zu(k) - inflow_damping_height ) / &
	1275	inflow_damping_width
[767]	1276	ELSE
[1340]	1277	inflow_damping_factor(k) = 0.0_wp
[767]	1278	ENDIF
[151]	1279
[767]	1280	ENDDO
[151]	1281
[147]	1282	ENDIF
	1283
[152]	1284	!
[2696]	1285	!-- Inside buildings set velocities back to zero
[1788]	1286	IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND. &
[359]	1287	topography /= 'flat' ) THEN
	1288	!
[2696]	1289	!-- Inside buildings set velocities back to zero.
	1290	!-- Other scalars (pt, q, s, p, sa, ...) are ignored at present,
[359]	1291	!-- maybe revise later.
[1001]	1292	DO i = nxlg, nxrg
	1293	DO j = nysg, nyng
[2232]	1294	DO k = nzb, nzt
	1295	u(k,j,i) = MERGE( u(k,j,i), 0.0_wp, &
	1296	BTEST( wall_flags_0(k,j,i), 1 ) )
	1297	v(k,j,i) = MERGE( v(k,j,i), 0.0_wp, &
	1298	BTEST( wall_flags_0(k,j,i), 2 ) )
	1299	w(k,j,i) = MERGE( w(k,j,i), 0.0_wp, &
	1300	BTEST( wall_flags_0(k,j,i), 3 ) )
	1301	ENDDO
[359]	1302	ENDDO
[1001]	1303	ENDDO
[359]	1304
	1305	ENDIF
	1306
	1307	!
[1]	1308	!-- Calculate initial temperature field and other constants used in case
	1309	!-- of a sloping surface
	1310	IF ( sloping_surface ) CALL init_slope
	1311
	1312	!
	1313	!-- Initialize new time levels (only done in order to set boundary values
	1314	!-- including ghost points)
[2696]	1315	pt_p = pt; u_p = u; v_p = v; w_p = w
[1960]	1316	IF ( humidity ) THEN
[1053]	1317	q_p = q
	1318	ENDIF
[1960]	1319	IF ( passive_scalar ) s_p = s
[181]	1320	!
	1321	!-- Allthough tendency arrays are set in prognostic_equations, they have
	1322	!-- have to be predefined here because they are used (but multiplied with 0)
	1323	!-- there before they are set.
[2696]	1324	tpt_m = 0.0_wp; tu_m = 0.0_wp; tv_m = 0.0_wp; tw_m = 0.0_wp
[1960]	1325	IF ( humidity ) THEN
[1340]	1326	tq_m = 0.0_wp
[1053]	1327	ENDIF
[1960]	1328	IF ( passive_scalar ) ts_m = 0.0_wp
[181]	1329
[3987]	1330	IF ( debug_output ) CALL debug_message( 'initializing in case of restart / cyclic_fill', 'end' )
[1384]	1331
[1]	1332	ELSE
	1333	!
	1334	!-- Actually this part of the programm should not be reached
[254]	1335	message_string = 'unknown initializing problem'
	1336	CALL message( 'init_3d_model', 'PA0193', 1, 2, 0, 6, 0 )
[1]	1337	ENDIF
	1338
[151]	1339
	1340	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
[1]	1341	!
[151]	1342	!-- Initialize old timelevels needed for radiation boundary conditions
[3182]	1343	IF ( bc_radiation_l ) THEN
[151]	1344	u_m_l(:,:,:) = u(:,:,1:2)
	1345	v_m_l(:,:,:) = v(:,:,0:1)
	1346	w_m_l(:,:,:) = w(:,:,0:1)
	1347	ENDIF
[3182]	1348	IF ( bc_radiation_r ) THEN
[151]	1349	u_m_r(:,:,:) = u(:,:,nx-1:nx)
	1350	v_m_r(:,:,:) = v(:,:,nx-1:nx)
	1351	w_m_r(:,:,:) = w(:,:,nx-1:nx)
	1352	ENDIF
[3182]	1353	IF ( bc_radiation_s ) THEN
[151]	1354	u_m_s(:,:,:) = u(:,0:1,:)
	1355	v_m_s(:,:,:) = v(:,1:2,:)
	1356	w_m_s(:,:,:) = w(:,0:1,:)
	1357	ENDIF
[3182]	1358	IF ( bc_radiation_n ) THEN
[151]	1359	u_m_n(:,:,:) = u(:,ny-1:ny,:)
	1360	v_m_n(:,:,:) = v(:,ny-1:ny,:)
	1361	w_m_n(:,:,:) = w(:,ny-1:ny,:)
	1362	ENDIF
[667]	1363
[151]	1364	ENDIF
[680]	1365
[667]	1366	!
	1367	!-- Calculate the initial volume flow at the right and north boundary
[709]	1368	IF ( conserve_volume_flow ) THEN
[151]	1369
[767]	1370	IF ( use_prescribed_profile_data ) THEN
[667]	1371
[1340]	1372	volume_flow_initial_l = 0.0_wp
	1373	volume_flow_area_l = 0.0_wp
[732]	1374
[667]	1375	IF ( nxr == nx ) THEN
	1376	DO j = nys, nyn
[2232]	1377	DO k = nzb+1, nzt
[1788]	1378	volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
[2232]	1379	u_init(k) * dzw(k) &
	1380	* MERGE( 1.0_wp, 0.0_wp, &
	1381	BTEST( wall_flags_0(k,j,nxr), 1 )&
	1382	)
	1383
	1384	volume_flow_area_l(1) = volume_flow_area_l(1) + dzw(k) &
	1385	* MERGE( 1.0_wp, 0.0_wp, &
	1386	BTEST( wall_flags_0(k,j,nxr), 1 )&
	1387	)
[767]	1388	ENDDO
	1389	ENDDO
	1390	ENDIF
	1391
	1392	IF ( nyn == ny ) THEN
	1393	DO i = nxl, nxr
[2232]	1394	DO k = nzb+1, nzt
	1395	volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
	1396	v_init(k) * dzw(k) &
	1397	* MERGE( 1.0_wp, 0.0_wp, &
	1398	BTEST( wall_flags_0(k,nyn,i), 2 )&
	1399	)
	1400	volume_flow_area_l(2) = volume_flow_area_l(2) + dzw(k) &
	1401	* MERGE( 1.0_wp, 0.0_wp, &
	1402	BTEST( wall_flags_0(k,nyn,i), 2 )&
	1403	)
[767]	1404	ENDDO
	1405	ENDDO
	1406	ENDIF
	1407
	1408	#if defined( __parallel )
	1409	CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&
	1410	2, MPI_REAL, MPI_SUM, comm2d, ierr )
	1411	CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1), &
	1412	2, MPI_REAL, MPI_SUM, comm2d, ierr )
	1413
	1414	#else
	1415	volume_flow_initial = volume_flow_initial_l
	1416	volume_flow_area = volume_flow_area_l
	1417	#endif
	1418
	1419	ELSEIF ( TRIM( initializing_actions ) == 'cyclic_fill' ) THEN
	1420
[1340]	1421	volume_flow_initial_l = 0.0_wp
	1422	volume_flow_area_l = 0.0_wp
[767]	1423
	1424	IF ( nxr == nx ) THEN
	1425	DO j = nys, nyn
[2232]	1426	DO k = nzb+1, nzt
[1788]	1427	volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
[2232]	1428	hom_sum(k,1,0) * dzw(k) &
	1429	* MERGE( 1.0_wp, 0.0_wp, &
	1430	BTEST( wall_flags_0(k,j,nx), 1 ) &
	1431	)
	1432	volume_flow_area_l(1) = volume_flow_area_l(1) + dzw(k) &
	1433	* MERGE( 1.0_wp, 0.0_wp, &
	1434	BTEST( wall_flags_0(k,j,nx), 1 ) &
	1435	)
[667]	1436	ENDDO
	1437	ENDDO
	1438	ENDIF
	1439
	1440	IF ( nyn == ny ) THEN
	1441	DO i = nxl, nxr
[2232]	1442	DO k = nzb+1, nzt
[1788]	1443	volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
[2232]	1444	hom_sum(k,2,0) * dzw(k) &
	1445	* MERGE( 1.0_wp, 0.0_wp, &
	1446	BTEST( wall_flags_0(k,ny,i), 2 ) &
	1447	)
	1448	volume_flow_area_l(2) = volume_flow_area_l(2) + dzw(k) &
	1449	* MERGE( 1.0_wp, 0.0_wp, &
	1450	BTEST( wall_flags_0(k,ny,i), 2 ) &
	1451	)
[667]	1452	ENDDO
	1453	ENDDO
	1454	ENDIF
	1455
[732]	1456	#if defined( __parallel )
	1457	CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&
	1458	2, MPI_REAL, MPI_SUM, comm2d, ierr )
	1459	CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1), &
	1460	2, MPI_REAL, MPI_SUM, comm2d, ierr )
	1461
	1462	#else
	1463	volume_flow_initial = volume_flow_initial_l
	1464	volume_flow_area = volume_flow_area_l
	1465	#endif
	1466
[667]	1467	ELSEIF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
	1468
[1340]	1469	volume_flow_initial_l = 0.0_wp
	1470	volume_flow_area_l = 0.0_wp
[732]	1471
[667]	1472	IF ( nxr == nx ) THEN
	1473	DO j = nys, nyn
[2232]	1474	DO k = nzb+1, nzt
	1475	volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
	1476	u(k,j,nx) * dzw(k) &
	1477	* MERGE( 1.0_wp, 0.0_wp, &
	1478	BTEST( wall_flags_0(k,j,nx), 1 ) &
	1479	)
	1480	volume_flow_area_l(1) = volume_flow_area_l(1) + dzw(k) &
	1481	* MERGE( 1.0_wp, 0.0_wp, &
	1482	BTEST( wall_flags_0(k,j,nx), 1 ) &
	1483	)
[667]	1484	ENDDO
	1485	ENDDO
	1486	ENDIF
	1487
	1488	IF ( nyn == ny ) THEN
	1489	DO i = nxl, nxr
[2232]	1490	DO k = nzb+1, nzt
[1788]	1491	volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
[2232]	1492	v(k,ny,i) * dzw(k) &
	1493	* MERGE( 1.0_wp, 0.0_wp, &
	1494	BTEST( wall_flags_0(k,ny,i), 2 ) &
	1495	)
	1496	volume_flow_area_l(2) = volume_flow_area_l(2) + dzw(k) &
	1497	* MERGE( 1.0_wp, 0.0_wp, &
	1498	BTEST( wall_flags_0(k,ny,i), 2 ) &
	1499	)
[667]	1500	ENDDO
	1501	ENDDO
	1502	ENDIF
	1503
	1504	#if defined( __parallel )
[732]	1505	CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&
	1506	2, MPI_REAL, MPI_SUM, comm2d, ierr )
	1507	CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1), &
	1508	2, MPI_REAL, MPI_SUM, comm2d, ierr )
[667]	1509
	1510	#else
[732]	1511	volume_flow_initial = volume_flow_initial_l
	1512	volume_flow_area = volume_flow_area_l
[667]	1513	#endif
	1514
[732]	1515	ENDIF
	1516
[151]	1517	!
[709]	1518	!-- In case of 'bulk_velocity' mode, volume_flow_initial is calculated
	1519	!-- from u\|v_bulk instead
[680]	1520	IF ( TRIM( conserve_volume_flow_mode ) == 'bulk_velocity' ) THEN
	1521	volume_flow_initial(1) = u_bulk * volume_flow_area(1)
	1522	volume_flow_initial(2) = v_bulk * volume_flow_area(2)
	1523	ENDIF
[667]	1524
[680]	1525	ENDIF
[2232]	1526	!
[4150]	1527	!-- In the following, surface properties can be further initialized with
	1528	!-- input from static driver file.
	1529	!-- At the moment this affects only default surfaces. For example,
	1530	!-- roughness length or sensible / latent heat fluxes can be initialized
	1531	!-- heterogeneously for default surfaces. Therefore, a generic routine
	1532	!-- from netcdf_data_input_mod is called to read a 2D array.
	1533	IF ( input_pids_static ) THEN
	1534	!
[4151]	1535	!-- Allocate memory for possible static input
	1536	ALLOCATE( tmp_2d%var(nys:nyn,nxl:nxr) )
	1537	tmp_2d%var = 0.0_wp
	1538	!
[4150]	1539	!-- Open the static input file
[4151]	1540	#if defined( __netcdf )
[4150]	1541	CALL open_read_file( TRIM( input_file_static ) // &
	1542	TRIM( coupling_char ), &
	1543	pids_static_id )
	1544
	1545	CALL inquire_num_variables( pids_static_id, num_vars_pids_static )
	1546	!
	1547	!-- Allocate memory to store variable names and read them
	1548	ALLOCATE( vars_pids_static(1:num_vars_pids_static) )
	1549	CALL inquire_variable_names( pids_static_id, vars_pids_static )
	1550	!
	1551	!-- Input roughness length.
	1552	IF ( check_existence( vars_pids_static, 'z0' ) ) THEN
	1553	!
	1554	!-- Read _FillValue attribute
	1555	CALL get_attribute( pids_static_id, char_fill, tmp_2d%fill, &
	1556	.FALSE., 'z0' )
	1557	!
	1558	!-- Read variable
	1559	CALL get_variable( pids_static_id, 'z0', tmp_2d%var, &
	1560	nxl, nxr, nys, nyn )
	1561	!
	1562	!-- Initialize roughness length. Note, z0 will be only initialized at
	1563	!-- default-type surfaces. At natural or urban z0 is implicitly
	1564	!-- initialized bythe respective parameter lists.
	1565	!-- Initialize horizontal surface elements.
	1566	CALL init_single_surface_properties( surf_def_h(0)%z0, &
	1567	tmp_2d%var, &
	1568	surf_def_h(0)%ns, &
	1569	tmp_2d%fill, &
	1570	surf_def_h(0)%i, &
	1571	surf_def_h(0)%j )
	1572	!
	1573	!-- Initialize roughness also at vertical surface elements.
	1574	!-- Note, the actual 2D input arrays are only defined on the
	1575	!-- subdomain. Therefore, pass the index arrays with their respective
	1576	!-- offset values.
	1577	DO l = 0, 3
	1578	CALL init_single_surface_properties( &
	1579	surf_def_v(l)%z0, &
	1580	tmp_2d%var, &
	1581	surf_def_v(l)%ns, &
	1582	tmp_2d%fill, &
	1583	surf_def_v(l)%i + surf_def_v(l)%ioff, &
	1584	surf_def_v(l)%j + surf_def_v(l)%joff )
	1585	ENDDO
	1586
	1587	ENDIF
	1588	!
	1589	!-- Additional variables, e.g. shf, qsws, etc, can be initialized the
	1590	!-- same way.
	1591
	1592	!
	1593	!-- Finally, close the input file.
	1594	CALL close_input_file( pids_static_id )
[4151]	1595	#endif
	1596	DEALLOCATE( tmp_2d%var )
[4150]	1597	ENDIF
	1598	!
[2618]	1599	!-- Finally, if random_heatflux is set, disturb shf at horizontal
	1600	!-- surfaces. Actually, this should be done in surface_mod, where all other
	1601	!-- initializations of surface quantities are done. However, this
	1602	!-- would create a ring dependency, hence, it is done here. Maybe delete
	1603	!-- disturb_heatflux and tranfer the respective code directly into the
	1604	!-- initialization in surface_mod.
[2232]	1605	IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. &
	1606	TRIM( initializing_actions ) /= 'cyclic_fill' ) THEN
[2618]	1607
[2232]	1608	IF ( use_surface_fluxes .AND. constant_heatflux .AND. &
	1609	random_heatflux ) THEN
	1610	IF ( surf_def_h(0)%ns >= 1 ) CALL disturb_heatflux( surf_def_h(0) )
	1611	IF ( surf_lsm_h%ns >= 1 ) CALL disturb_heatflux( surf_lsm_h )
	1612	IF ( surf_usm_h%ns >= 1 ) CALL disturb_heatflux( surf_usm_h )
	1613	ENDIF
	1614	ENDIF
[680]	1615
[787]	1616	!
[3747]	1617	!-- Compute total sum of grid points and the mean surface level height for each
	1618	!-- statistic region. These are mainly used for horizontal averaging of
	1619	!-- turbulence statistics.
[2696]	1620	!-- ngp_2dh: number of grid points of a horizontal cross section through the
[3747]	1621	!-- respective statistic region
	1622	!-- ngp_3d: number of grid points of the respective statistic region
[2696]	1623	ngp_2dh_outer_l = 0
	1624	ngp_2dh_outer = 0
	1625	ngp_2dh_s_inner_l = 0
	1626	ngp_2dh_s_inner = 0
	1627	ngp_2dh_l = 0
	1628	ngp_2dh = 0
	1629	ngp_3d_inner_l = 0.0_wp
	1630	ngp_3d_inner = 0
	1631	ngp_3d = 0
	1632	ngp_sums = ( nz + 2 ) * ( pr_palm + max_pr_user )
	1633
	1634	mean_surface_level_height = 0.0_wp
	1635	mean_surface_level_height_l = 0.0_wp
	1636	!
	1637	!-- To do: New concept for these non-topography grid points!
	1638	DO sr = 0, statistic_regions
	1639	DO i = nxl, nxr
	1640	DO j = nys, nyn
	1641	IF ( rmask(j,i,sr) == 1.0_wp ) THEN
	1642	!
	1643	!-- All xy-grid points
	1644	ngp_2dh_l(sr) = ngp_2dh_l(sr) + 1
	1645	!
	1646	!-- Determine mean surface-level height. In case of downward-
	1647	!-- facing walls are present, more than one surface level exist.
	1648	!-- In this case, use the lowest surface-level height.
	1649	IF ( surf_def_h(0)%start_index(j,i) <= &
	1650	surf_def_h(0)%end_index(j,i) ) THEN
	1651	m = surf_def_h(0)%start_index(j,i)
	1652	k = surf_def_h(0)%k(m)
	1653	mean_surface_level_height_l(sr) = &
	1654	mean_surface_level_height_l(sr) + zw(k-1)
	1655	ENDIF
	1656	IF ( surf_lsm_h%start_index(j,i) <= &
	1657	surf_lsm_h%end_index(j,i) ) THEN
	1658	m = surf_lsm_h%start_index(j,i)
	1659	k = surf_lsm_h%k(m)
	1660	mean_surface_level_height_l(sr) = &
	1661	mean_surface_level_height_l(sr) + zw(k-1)
	1662	ENDIF
	1663	IF ( surf_usm_h%start_index(j,i) <= &
	1664	surf_usm_h%end_index(j,i) ) THEN
	1665	m = surf_usm_h%start_index(j,i)
	1666	k = surf_usm_h%k(m)
	1667	mean_surface_level_height_l(sr) = &
	1668	mean_surface_level_height_l(sr) + zw(k-1)
	1669	ENDIF
	1670
	1671	k_surf = k - 1
	1672
	1673	DO k = nzb, nzt+1
	1674	!
	1675	!-- xy-grid points above topography
	1676	ngp_2dh_outer_l(k,sr) = ngp_2dh_outer_l(k,sr) + &
	1677	MERGE( 1, 0, BTEST( wall_flags_0(k,j,i), 24 ) )
	1678
	1679	ngp_2dh_s_inner_l(k,sr) = ngp_2dh_s_inner_l(k,sr) + &
	1680	MERGE( 1, 0, BTEST( wall_flags_0(k,j,i), 22 ) )
	1681
	1682	ENDDO
	1683	!
	1684	!-- All grid points of the total domain above topography
	1685	ngp_3d_inner_l(sr) = ngp_3d_inner_l(sr) + ( nz - k_surf + 2 )
	1686
	1687
	1688
	1689	ENDIF
	1690	ENDDO
	1691	ENDDO
	1692	ENDDO
[3747]	1693
[2696]	1694	sr = statistic_regions + 1
	1695	#if defined( __parallel )
	1696	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
	1697	CALL MPI_ALLREDUCE( ngp_2dh_l(0), ngp_2dh(0), sr, MPI_INTEGER, MPI_SUM, &
	1698	comm2d, ierr )
	1699	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
	1700	CALL MPI_ALLREDUCE( ngp_2dh_outer_l(0,0), ngp_2dh_outer(0,0), (nz+2)*sr, &
	1701	MPI_INTEGER, MPI_SUM, comm2d, ierr )
	1702	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
	1703	CALL MPI_ALLREDUCE( ngp_2dh_s_inner_l(0,0), ngp_2dh_s_inner(0,0), &
	1704	(nz+2)*sr, MPI_INTEGER, MPI_SUM, comm2d, ierr )
	1705	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
	1706	CALL MPI_ALLREDUCE( ngp_3d_inner_l(0), ngp_3d_inner_tmp(0), sr, MPI_REAL, &
	1707	MPI_SUM, comm2d, ierr )
	1708	ngp_3d_inner = INT( ngp_3d_inner_tmp, KIND = SELECTED_INT_KIND( 18 ) )
	1709	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
	1710	CALL MPI_ALLREDUCE( mean_surface_level_height_l(0), &
	1711	mean_surface_level_height(0), sr, MPI_REAL, &
	1712	MPI_SUM, comm2d, ierr )
	1713	mean_surface_level_height = mean_surface_level_height / REAL( ngp_2dh )
	1714	#else
	1715	ngp_2dh = ngp_2dh_l
	1716	ngp_2dh_outer = ngp_2dh_outer_l
	1717	ngp_2dh_s_inner = ngp_2dh_s_inner_l
	1718	ngp_3d_inner = INT( ngp_3d_inner_l, KIND = SELECTED_INT_KIND( 18 ) )
	1719	mean_surface_level_height = mean_surface_level_height_l / REAL( ngp_2dh_l )
	1720	#endif
	1721
	1722	ngp_3d = INT ( ngp_2dh, KIND = SELECTED_INT_KIND( 18 ) ) * &
	1723	INT ( (nz + 2 ), KIND = SELECTED_INT_KIND( 18 ) )
	1724
	1725	!
	1726	!-- Set a lower limit of 1 in order to avoid zero divisions in flow_statistics,
	1727	!-- buoyancy, etc. A zero value will occur for cases where all grid points of
	1728	!-- the respective subdomain lie below the surface topography
	1729	ngp_2dh_outer = MAX( 1, ngp_2dh_outer(:,:) )
	1730	ngp_3d_inner = MAX( INT(1, KIND = SELECTED_INT_KIND( 18 )), &
	1731	ngp_3d_inner(:) )
	1732	ngp_2dh_s_inner = MAX( 1, ngp_2dh_s_inner(:,:) )
	1733
	1734	DEALLOCATE( mean_surface_level_height_l, ngp_2dh_l, ngp_2dh_outer_l, &
	1735	ngp_3d_inner_l, ngp_3d_inner_tmp )
	1736	!
[2232]	1737	!-- Initialize surface forcing corresponding to large-scale forcing. Therein,
	1738	!-- initialize heat-fluxes, etc. via datatype. Revise it later!
	1739	IF ( large_scale_forcing .AND. lsf_surf ) THEN
	1740	IF ( use_surface_fluxes .AND. constant_heatflux ) THEN
	1741	CALL ls_forcing_surf ( simulated_time )
	1742	ENDIF
	1743	ENDIF
	1744	!
[3347]	1745	!-- Initializae 3D offline nesting in COSMO model and read data from
	1746	!-- external NetCDF file.
	1747	IF ( nesting_offline ) CALL nesting_offl_init
	1748	!
[787]	1749	!-- Initialize quantities for special advections schemes
	1750	CALL init_advec
[680]	1751
[667]	1752	!
[680]	1753	!-- Impose random perturbation on the horizontal velocity field and then
	1754	!-- remove the divergences from the velocity field at the initial stage
[1788]	1755	IF ( create_disturbances .AND. disturbance_energy_limit /= 0.0_wp .AND. &
	1756	TRIM( initializing_actions ) /= 'read_restart_data' .AND. &
[680]	1757	TRIM( initializing_actions ) /= 'cyclic_fill' ) THEN
	1758
[3987]	1759	IF ( debug_output ) CALL debug_message( 'creating disturbances + applying pressure solver', 'start' )
[3849]	1760	!
	1761	!-- Needed for both disturb_field and pres
	1762	!$ACC DATA &
	1763	!$ACC CREATE(tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
	1764	!$ACC COPY(u(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
	1765	!$ACC COPY(v(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
	1766
[2232]	1767	CALL disturb_field( 'u', tend, u )
	1768	CALL disturb_field( 'v', tend, v )
[1384]	1769
[3849]	1770	!$ACC DATA &
	1771	!$ACC CREATE(d(nzb+1:nzt,nys:nyn,nxl:nxr)) &
	1772	!$ACC COPY(w(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
	1773	!$ACC COPY(p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
	1774	!$ACC COPYIN(rho_air(nzb:nzt+1), rho_air_zw(nzb:nzt+1)) &
	1775	!$ACC COPYIN(ddzu(1:nzt+1), ddzw(1:nzt+1)) &
	1776	!$ACC COPYIN(wall_flags_0(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
	1777	!$ACC COPYIN(bc_h(0:1)) &
	1778	!$ACC COPYIN(bc_h(0)%i(1:bc_h(0)%ns)) &
	1779	!$ACC COPYIN(bc_h(0)%j(1:bc_h(0)%ns)) &
	1780	!$ACC COPYIN(bc_h(0)%k(1:bc_h(0)%ns)) &
	1781	!$ACC COPYIN(bc_h(1)%i(1:bc_h(1)%ns)) &
	1782	!$ACC COPYIN(bc_h(1)%j(1:bc_h(1)%ns)) &
	1783	!$ACC COPYIN(bc_h(1)%k(1:bc_h(1)%ns))
	1784
[680]	1785	n_sor = nsor_ini
	1786	CALL pres
	1787	n_sor = nsor
[1384]	1788
[3849]	1789	!$ACC END DATA
	1790	!$ACC END DATA
	1791
[3987]	1792	IF ( debug_output ) CALL debug_message( 'creating disturbances + applying pressure solver', 'end' )
	1793
[680]	1794	ENDIF
	1795
[3294]	1796	IF ( .NOT. ocean_mode ) THEN
[3274]	1797
	1798	ALLOCATE( hyp(nzb:nzt+1) )
	1799	ALLOCATE( d_exner(nzb:nzt+1) )
	1800	ALLOCATE( exner(nzb:nzt+1) )
	1801	ALLOCATE( hyrho(nzb:nzt+1) )
[1849]	1802	!
[3274]	1803	!-- Check temperature in case of too large domain height
	1804	DO k = nzb, nzt+1
	1805	IF ( ( pt_surface * exner_function(surface_pressure * 100.0_wp) - g/c_p * zu(k) ) < 0.0_wp ) THEN
	1806	WRITE( message_string, * ) 'absolute temperature < 0.0 at zu(', k, &
	1807	') = ', zu(k)
[3685]	1808	CALL message( 'init_3d_model', 'PA0142', 1, 2, 0, 6, 0 )
[3274]	1809	ENDIF
	1810	ENDDO
	1811
	1812	!
	1813	!-- Calculate vertical profile of the hydrostatic pressure (hyp)
	1814	hyp = barometric_formula(zu, pt_surface * exner_function(surface_pressure * 100.0_wp), surface_pressure * 100.0_wp)
	1815	d_exner = exner_function_invers(hyp)
	1816	exner = 1.0_wp / exner_function_invers(hyp)
	1817	hyrho = ideal_gas_law_rho_pt(hyp, pt_init)
	1818	!
	1819	!-- Compute reference density
	1820	rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp, pt_surface * exner_function(surface_pressure * 100.0_wp))
	1821
[96]	1822	ENDIF
[1]	1823
	1824	!
	1825	!-- If required, initialize particles
[3159]	1826	IF ( agents_active ) CALL mas_init
	1827	!
[3937]	1828	!-- Initialization of synthetic turbulence generator
	1829	IF ( use_syn_turb_gen ) CALL stg_init
[2696]	1830	!
[3685]	1831	!-- Initializing actions for all other modules
	1832	CALL module_interface_init
[2696]	1833	!
[3685]	1834	!-- Initialize surface layer (done after LSM as roughness length are required
	1835	!-- for initialization
	1836	IF ( constant_flux_layer ) CALL init_surface_layer_fluxes
[2977]	1837	!
[3421]	1838	!-- Initialize surface data output
[3685]	1839	IF ( surface_output ) CALL surface_data_output_init
[3472]	1840	!
[673]	1841	!-- Initialize the ws-scheme.
[3448]	1842	IF ( ws_scheme_sca .OR. ws_scheme_mom ) CALL ws_init
[3711]	1843	!
	1844	!-- Perform post-initializing checks for all other modules
	1845	CALL module_interface_init_checks
[1]	1846
	1847	!
[709]	1848	!-- Setting weighting factors for calculation of perturbation pressure
[1762]	1849	!-- and turbulent quantities from the RK substeps
[709]	1850	IF ( TRIM(timestep_scheme) == 'runge-kutta-3' ) THEN ! for RK3-method
	1851
[1322]	1852	weight_substep(1) = 1._wp/6._wp
	1853	weight_substep(2) = 3._wp/10._wp
	1854	weight_substep(3) = 8._wp/15._wp
[709]	1855
[1322]	1856	weight_pres(1) = 1._wp/3._wp
	1857	weight_pres(2) = 5._wp/12._wp
	1858	weight_pres(3) = 1._wp/4._wp
[709]	1859
	1860	ELSEIF ( TRIM(timestep_scheme) == 'runge-kutta-2' ) THEN ! for RK2-method
	1861
[1322]	1862	weight_substep(1) = 1._wp/2._wp
	1863	weight_substep(2) = 1._wp/2._wp
[673]	1864
[1322]	1865	weight_pres(1) = 1._wp/2._wp
	1866	weight_pres(2) = 1._wp/2._wp
[709]	1867
[1001]	1868	ELSE ! for Euler-method
[709]	1869
[1340]	1870	weight_substep(1) = 1.0_wp
	1871	weight_pres(1) = 1.0_wp
[709]	1872
[673]	1873	ENDIF
	1874
	1875	!
[1]	1876	!-- Initialize Rayleigh damping factors
[1340]	1877	rdf = 0.0_wp
	1878	rdf_sc = 0.0_wp
	1879	IF ( rayleigh_damping_factor /= 0.0_wp ) THEN
[3294]	1880
	1881	IF ( .NOT. ocean_mode ) THEN
[108]	1882	DO k = nzb+1, nzt
	1883	IF ( zu(k) >= rayleigh_damping_height ) THEN
[1788]	1884	rdf(k) = rayleigh_damping_factor * &
[1340]	1885	( SIN( pi * 0.5_wp * ( zu(k) - rayleigh_damping_height ) &
[1788]	1886	/ ( zu(nzt) - rayleigh_damping_height ) ) &
[1]	1887	)**2
[108]	1888	ENDIF
	1889	ENDDO
	1890	ELSE
[3294]	1891	!
	1892	!-- In ocean mode, rayleigh damping is applied in the lower part of the
	1893	!-- model domain
[108]	1894	DO k = nzt, nzb+1, -1
	1895	IF ( zu(k) <= rayleigh_damping_height ) THEN
[1788]	1896	rdf(k) = rayleigh_damping_factor * &
[1340]	1897	( SIN( pi * 0.5_wp * ( rayleigh_damping_height - zu(k) ) &
[1788]	1898	/ ( rayleigh_damping_height - zu(nzb+1) ) ) &
[108]	1899	)**2
	1900	ENDIF
	1901	ENDDO
	1902	ENDIF
[3294]	1903
[1]	1904	ENDIF
[785]	1905	IF ( scalar_rayleigh_damping ) rdf_sc = rdf
[1]	1906
	1907	!
[240]	1908	!-- Initialize the starting level and the vertical smoothing factor used for
	1909	!-- the external pressure gradient
[1340]	1910	dp_smooth_factor = 1.0_wp
[240]	1911	IF ( dp_external ) THEN
	1912	!
	1913	!-- Set the starting level dp_level_ind_b only if it has not been set before
	1914	!-- (e.g. in init_grid).
	1915	IF ( dp_level_ind_b == 0 ) THEN
	1916	ind_array = MINLOC( ABS( dp_level_b - zu ) )
	1917	dp_level_ind_b = ind_array(1) - 1 + nzb
	1918	! MINLOC uses lower array bound 1
	1919	ENDIF
	1920	IF ( dp_smooth ) THEN
[1340]	1921	dp_smooth_factor(:dp_level_ind_b) = 0.0_wp
[240]	1922	DO k = dp_level_ind_b+1, nzt
[1340]	1923	dp_smooth_factor(k) = 0.5_wp * ( 1.0_wp + SIN( pi * &
	1924	( REAL( k - dp_level_ind_b, KIND=wp ) / &
	1925	REAL( nzt - dp_level_ind_b, KIND=wp ) - 0.5_wp ) ) )
[240]	1926	ENDDO
	1927	ENDIF
	1928	ENDIF
	1929
	1930	!
[978]	1931	!-- Initialize damping zone for the potential temperature in case of
	1932	!-- non-cyclic lateral boundaries. The damping zone has the maximum value
	1933	!-- at the inflow boundary and decreases to zero at pt_damping_width.
[1340]	1934	ptdf_x = 0.0_wp
	1935	ptdf_y = 0.0_wp
[1159]	1936	IF ( bc_lr_dirrad ) THEN
[996]	1937	DO i = nxl, nxr
[978]	1938	IF ( ( i * dx ) < pt_damping_width ) THEN
[1340]	1939	ptdf_x(i) = pt_damping_factor * ( SIN( pi * 0.5_wp * &
	1940	REAL( pt_damping_width - i * dx, KIND=wp ) / ( &
[1788]	1941	REAL( pt_damping_width, KIND=wp ) ) ) )**2
[73]	1942	ENDIF
	1943	ENDDO
[1159]	1944	ELSEIF ( bc_lr_raddir ) THEN
[996]	1945	DO i = nxl, nxr
[978]	1946	IF ( ( i * dx ) > ( nx * dx - pt_damping_width ) ) THEN
[1322]	1947	ptdf_x(i) = pt_damping_factor * &
[1340]	1948	SIN( pi * 0.5_wp * &
	1949	( ( i - nx ) * dx + pt_damping_width ) / &
	1950	REAL( pt_damping_width, KIND=wp ) )**2
[73]	1951	ENDIF
[978]	1952	ENDDO
[1159]	1953	ELSEIF ( bc_ns_dirrad ) THEN
[996]	1954	DO j = nys, nyn
[978]	1955	IF ( ( j * dy ) > ( ny * dy - pt_damping_width ) ) THEN
[1322]	1956	ptdf_y(j) = pt_damping_factor * &
[1340]	1957	SIN( pi * 0.5_wp * &
	1958	( ( j - ny ) * dy + pt_damping_width ) / &
	1959	REAL( pt_damping_width, KIND=wp ) )**2
[1]	1960	ENDIF
[978]	1961	ENDDO
[1159]	1962	ELSEIF ( bc_ns_raddir ) THEN
[996]	1963	DO j = nys, nyn
[978]	1964	IF ( ( j * dy ) < pt_damping_width ) THEN
[1322]	1965	ptdf_y(j) = pt_damping_factor * &
[1340]	1966	SIN( pi * 0.5_wp * &
	1967	( pt_damping_width - j * dy ) / &
	1968	REAL( pt_damping_width, KIND=wp ) )**2
[1]	1969	ENDIF
[73]	1970	ENDDO
[1]	1971	ENDIF
[51]	1972
[1]	1973	!
	1974	!-- Input binary data file is not needed anymore. This line must be placed
	1975	!-- after call of user_init!
	1976	CALL close_file( 13 )
[2934]	1977	!
	1978	!-- In case of nesting, put an barrier to assure that all parent and child
	1979	!-- domains finished initialization.
	1980	#if defined( __parallel )
	1981	IF ( nested_run ) CALL MPI_BARRIER( MPI_COMM_WORLD, ierr )
	1982	#endif
[1]	1983
[2934]	1984
[3987]	1985	CALL location_message( 'model initialization', 'finished' )
[1]	1986
	1987	END SUBROUTINE init_3d_model

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

Context Navigation

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

Download in other formats: