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

Line
1	!> @file init_3d_model.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of the PALM model system.
4	!
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.
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	!
17	! Copyright 1997-2019 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! ------------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: init_3d_model.f90 4185 2019-08-23 13:49:38Z oliver.maas $
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
34	! Corrected "Former revisions" section
35	!
36	! 4168 2019-08-16 13:50:17Z suehring
37	! Replace function get_topography_top_index by topo_top_ind
38	!
39	! 4151 2019-08-09 08:24:30Z suehring
40	! Add netcdf directive around input calls (fix for last commit)
41	!
42	! 4150 2019-08-08 20:00:47Z suehring
43	! Input of additional surface variables independent on land- or urban-surface
44	! model
45	!
46	! 4131 2019-08-02 11:06:18Z monakurppa
47	! Allocate sums and sums_l to allow profile output for salsa variables.
48	!
49	! 4130 2019-08-01 13:04:13Z suehring
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
55	! Unused variables removed
56	!
57	! 4088 2019-07-11 13:57:56Z Giersch
58	! Pressure and density profile calculations revised using basic functions
59	!
60	! 4048 2019-06-21 21:00:21Z knoop
61	! Further modularization of particle code components
62	!
63	! 4017 2019-06-06 12:16:46Z schwenkel
64	! Convert most location messages to debug messages to reduce output in
65	! job logfile to a minimum
66	!
67	!
68	! unused variable removed
69	!
70	! 3937 2019-04-29 15:09:07Z suehring
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
76	! Fix problem with LOD = 2 initialization
77	!
78	! 3885 2019-04-11 11:29:34Z kanani
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
83	! Move initialization of rmask before initializing user_init_arrays
84	!
85	! 3711 2019-01-31 13:44:26Z knoop
86	! Introduced module_interface_init_checks for post-init checks in modules
87	!
88	! 3700 2019-01-26 17:03:42Z knoop
89	! Some interface calls moved to module_interface + cleanup
90	!
91	! 3648 2019-01-02 16:35:46Z suehring
92	! Rename subroutines for surface-data output
93	!
94	! Revision 1.1 1998/03/09 16:22:22 raasch
95	! Initial revision
96	!
97	!
98	! Description:
99	! ------------
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
106	!------------------------------------------------------------------------------!
107	SUBROUTINE init_3d_model
108
109
110	USE advec_ws
111
112	USE arrays_3d
113
114	USE basic_constants_and_equations_mod, &
115	ONLY: c_p, g, l_v, pi, exner_function, exner_function_invers, &
116	ideal_gas_law_rho, ideal_gas_law_rho_pt, barometric_formula
117
118	USE bulk_cloud_model_mod, &
119	ONLY: bulk_cloud_model
120
121	USE chem_modules, &
122	ONLY: max_pr_cs ! ToDo: this dependency needs to be removed cause it is ugly #new_dom
123
124	USE control_parameters
125
126	USE grid_variables, &
127	ONLY: dx, dy, ddx2_mg, ddy2_mg
128
129	USE indices
130
131	USE kinds
132
133	USE lsf_nudging_mod, &
134	ONLY: ls_forcing_surf
135
136	USE model_1d_mod, &
137	ONLY: init_1d_model, l1d, u1d, v1d
138
139	USE module_interface, &
140	ONLY: module_interface_init_arrays, &
141	module_interface_init, &
142	module_interface_init_checks
143
144	USE multi_agent_system_mod, &
145	ONLY: agents_active, mas_init
146
147	USE netcdf_interface, &
148	ONLY: dots_max
149
150	USE netcdf_data_input_mod, &
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
165	USE nesting_offl_mod, &
166	ONLY: nesting_offl_init
167
168	USE pegrid
169
170	#if defined( __parallel )
171	USE pmc_interface, &
172	ONLY: nested_run
173	#endif
174
175	USE random_function_mod
176
177	USE random_generator_parallel, &
178	ONLY: init_parallel_random_generator
179
180	USE read_restart_data_mod, &
181	ONLY: rrd_read_parts_of_global, rrd_local
182
183	USE statistics, &
184	ONLY: hom, hom_sum, mean_surface_level_height, pr_palm, rmask, &
185	statistic_regions, sums, sums_divnew_l, sums_divold_l, sums_l, &
186	sums_l_l, sums_wsts_bc_l, ts_value, &
187	weight_pres, weight_substep
188
189	USE synthetic_turbulence_generator_mod, &
190	ONLY: stg_init, use_syn_turb_gen
191
192	USE surface_layer_fluxes_mod, &
193	ONLY: init_surface_layer_fluxes
194
195	USE surface_mod, &
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, &
202	surf_usm_h
203
204	#if defined( _OPENACC )
205	USE surface_mod, &
206	ONLY : bc_h
207	#endif
208
209	USE surface_data_output_mod, &
210	ONLY: surface_data_output_init
211
212	USE transpose_indices
213
214	IMPLICIT NONE
215
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
237
238	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: ngp_2dh_l !< toal number of horizontal grid points in statistical region on subdomain
239
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
242
243
244
245	REAL(wp), DIMENSION(:), ALLOCATABLE :: init_l !< dummy array used for averaging 3D data to obtain inital profiles
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
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
253
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
257
258	TYPE(real_2d) :: tmp_2d !< temporary variable to input additional surface-data from static file
259
260	CALL location_message( 'model initialization', 'start' )
261
262	IF ( debug_output ) CALL debug_message( 'allocating arrays', 'start' )
263	!
264	!-- Allocate arrays
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), &
273	sums_divold_l(0:statistic_regions) )
274	ALLOCATE( dp_smooth_factor(nzb:nzt), rdf(nzb+1:nzt), rdf_sc(nzb+1:nzt) )
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), &
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), &
282	sums_l_l(nzb:nzt+1,0:statistic_regions,0:threads_per_task-1), &
283	sums_wsts_bc_l(nzb:nzt+1,0:statistic_regions) )
284	ALLOCATE( ts_value(dots_max,0:statistic_regions) )
285	ALLOCATE( ptdf_x(nxlg:nxrg), ptdf_y(nysg:nyng) )
286
287	ALLOCATE( d(nzb+1:nzt,nys:nyn,nxl:nxr), &
288	p(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
289	tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
290
291	ALLOCATE( pt_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
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), &
301	w_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
302	IF ( .NOT. neutral ) THEN
303	ALLOCATE( pt_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
304	ENDIF
305	!
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	!
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) )
319	ELSEIF ( psolver(1:9) == 'multigrid' ) THEN
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) )
323	ENDIF
324
325	!
326	!-- Array for storing constant coeffficients of the tridiagonal solver
327	IF ( psolver == 'poisfft' ) THEN
328	ALLOCATE( tri(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1,2) )
329	ALLOCATE( tric(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1) )
330	ENDIF
331
332	IF ( humidity ) THEN
333	!
334	!-- 3D-humidity
335	ALLOCATE( q_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
336	q_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
337	q_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
338	vpt_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
339	ENDIF
340
341	IF ( passive_scalar ) THEN
342
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) )
348
349	ENDIF
350
351	!
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	!
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	!
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
374	DO k = nzb, nzt+1
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))
380	ENDDO
381
382	DO k = nzb, nzt
383	rho_air_zw(k) = 0.5_wp * ( rho_air(k) + rho_air(k+1) )
384	ENDDO
385
386	rho_air_zw(nzt+1) = rho_air_zw(nzt) &
387	+ 2.0_wp * ( rho_air(nzt+1) - rho_air_zw(nzt) )
388
389	ELSE
390	DO k = nzb, nzt+1
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))
396	ENDDO
397
398	DO k = nzb, nzt
399	rho_air_zw(k) = 0.5_wp * ( rho_air(k) + rho_air(k+1) )
400	ENDDO
401
402	rho_air_zw(nzt+1) = rho_air_zw(nzt) &
403	+ 2.0_wp * ( rho_air(nzt+1) - rho_air_zw(nzt) )
404
405	ENDIF
406	!
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
429	heatflux_input_conversion(k) = 1.0_wp / c_p
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
439	heatflux_output_conversion(k) = c_p
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	!
513	!-- 1D-array for large scale subsidence velocity
514	IF ( .NOT. ALLOCATED( w_subs ) ) THEN
515	ALLOCATE ( w_subs(nzb:nzt+1) )
516	w_subs = 0.0_wp
517	ENDIF
518
519	!
520	!-- Arrays to store velocity data from t-dt and the phase speeds which
521	!-- are needed for radiation boundary conditions
522	IF ( bc_radiation_l ) THEN
523	ALLOCATE( u_m_l(nzb:nzt+1,nysg:nyng,1:2), &
524	v_m_l(nzb:nzt+1,nysg:nyng,0:1), &
525	w_m_l(nzb:nzt+1,nysg:nyng,0:1) )
526	ENDIF
527	IF ( bc_radiation_r ) THEN
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), &
530	w_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx) )
531	ENDIF
532	IF ( bc_radiation_l .OR. bc_radiation_r ) THEN
533	ALLOCATE( c_u(nzb:nzt+1,nysg:nyng), c_v(nzb:nzt+1,nysg:nyng), &
534	c_w(nzb:nzt+1,nysg:nyng) )
535	ENDIF
536	IF ( bc_radiation_s ) THEN
537	ALLOCATE( u_m_s(nzb:nzt+1,0:1,nxlg:nxrg), &
538	v_m_s(nzb:nzt+1,1:2,nxlg:nxrg), &
539	w_m_s(nzb:nzt+1,0:1,nxlg:nxrg) )
540	ENDIF
541	IF ( bc_radiation_n ) THEN
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), &
544	w_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg) )
545	ENDIF
546	IF ( bc_radiation_s .OR. bc_radiation_n ) THEN
547	ALLOCATE( c_u(nzb:nzt+1,nxlg:nxrg), c_v(nzb:nzt+1,nxlg:nxrg), &
548	c_w(nzb:nzt+1,nxlg:nxrg) )
549	ENDIF
550	IF ( bc_radiation_l .OR. bc_radiation_r .OR. bc_radiation_s .OR. &
551	bc_radiation_n ) THEN
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
555
556	!
557	!-- Initial assignment of the pointers
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
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
566
567	IF ( humidity ) THEN
568	q => q_1; q_p => q_2; tq_m => q_3
569	vpt => vpt_1
570	ENDIF
571
572	IF ( passive_scalar ) THEN
573	s => s_1; s_p => s_2; ts_m => s_3
574	ENDIF
575
576	!
577	!-- Initialize surface arrays
578	CALL init_surface_arrays
579	!
580	!-- Allocate arrays for other modules
581	CALL module_interface_init_arrays
582
583
584	!
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.
590	ALLOCATE( weight_substep(1:intermediate_timestep_count_max), &
591	weight_pres(1:intermediate_timestep_count_max) )
592	weight_substep = 1.0_wp
593	weight_pres = 1.0_wp
594	intermediate_timestep_count = 0 ! needed when simulated_time = 0.0
595
596	IF ( debug_output ) CALL debug_message( 'allocating arrays', 'end' )
597
598	!
599	!-- Initialize time series
600	ts_value = 0.0_wp
601
602	!
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
611
612	!
613	!-- Initialize model variables
614	IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. &
615	TRIM( initializing_actions ) /= 'cyclic_fill' ) THEN
616	!
617	!-- Initialization with provided input data derived from larger-scale model
618	IF ( INDEX( initializing_actions, 'inifor' ) /= 0 ) THEN
619	IF ( debug_output ) CALL debug_message( 'initializing with INIFOR', 'start' )
620	!
621	!-- Read initial 1D profiles or 3D data from NetCDF file, depending
622	!-- on the provided level-of-detail.
623	!-- At the moment, only u, v, w, pt and q are provided.
624	CALL netcdf_data_input_init_3d
625	!
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 )
641
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 )
649
650	ENDIF
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 )
660
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
710	!
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.
717	DO i = nxlg, nxrg
718	DO j = nysg, nyng
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 ) ) &
725	pt(:,j,i) = pt_init(:)
726	IF( humidity .AND. &
727	( init_3d%lod_q == 1 .OR. init_3d%lod_q == 2 ) ) &
728	q(:,j,i) = q_init(:)
729	ENDDO
730	ENDDO
731	!
732	!-- Set geostrophic wind components.
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
739	!
740	!-- Set bottom and top boundary condition for geostrophic wind
741	ug(nzt+1) = ug(nzt)
742	vg(nzt+1) = vg(nzt)
743	ug(nzb) = ug(nzb+1)
744	vg(nzb) = vg(nzb+1)
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 ) )
763	!
764	!-- Initialize surface variables, e.g. friction velocity, momentum
765	!-- fluxes, etc.
766	CALL init_surfaces
767
768	IF ( debug_output ) CALL debug_message( 'initializing with INIFOR', 'end' )
769	!
770	!-- Initialization via computed 1D-model profiles
771	ELSEIF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 ) THEN
772
773	IF ( debug_output ) CALL debug_message( 'initializing with 1D model profiles', 'start' )
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
780	DO i = nxlg, nxrg
781	DO j = nysg, nyng
782	pt(:,j,i) = pt_init
783	u(:,j,i) = u1d
784	v(:,j,i) = v1d
785	ENDDO
786	ENDDO
787
788	IF ( humidity ) THEN
789	DO i = nxlg, nxrg
790	DO j = nysg, nyng
791	q(:,j,i) = q_init
792	ENDDO
793	ENDDO
794	ENDIF
795
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
803	!
804	!-- Store initial profiles for output purposes etc.
805	IF ( .NOT. constant_diffusion ) THEN
806	hom(:,1,25,:) = SPREAD( l1d, 2, statistic_regions+1 )
807	ENDIF
808	!
809	!-- Set velocities back to zero
810	u = MERGE( u, 0.0_wp, BTEST( wall_flags_0, 1 ) )
811	v = MERGE( v, 0.0_wp, BTEST( wall_flags_0, 2 ) )
812	!
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
820	!
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)
826	ENDDO
827	ENDDO
828
829	ENDIF
830	!
831	!-- Initialize surface variables, e.g. friction velocity, momentum
832	!-- fluxes, etc.
833	CALL init_surfaces
834
835	IF ( debug_output ) CALL debug_message( 'initializing with 1D model profiles', 'end' )
836
837	ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 ) &
838	THEN
839
840	IF ( debug_output ) CALL debug_message( 'initializing with constant profiles', 'start' )
841
842	!
843	!-- Use constructed initial profiles (velocity constant with height,
844	!-- temperature profile with constant gradient)
845	DO i = nxlg, nxrg
846	DO j = nysg, nyng
847	pt(:,j,i) = pt_init
848	u(:,j,i) = u_init
849	v(:,j,i) = v_init
850	ENDDO
851	ENDDO
852	!
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	!
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
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
866	DO i = nxlg, nxrg
867	DO j = nysg, nyng
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
874	ENDDO
875	ENDDO
876	ENDIF
877
878	IF ( humidity ) THEN
879	DO i = nxlg, nxrg
880	DO j = nysg, nyng
881	q(:,j,i) = q_init
882	ENDDO
883	ENDDO
884	ENDIF
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
893
894	!
895	!-- Compute initial temperature field and other constants used in case
896	!-- of a sloping surface
897	IF ( sloping_surface ) CALL init_slope
898	!
899	!-- Initialize surface variables, e.g. friction velocity, momentum
900	!-- fluxes, etc.
901	CALL init_surfaces
902
903	IF ( debug_output ) CALL debug_message( 'initializing with constant profiles', 'end' )
904
905	ELSEIF ( INDEX(initializing_actions, 'by_user') /= 0 ) &
906	THEN
907
908	IF ( debug_output ) CALL debug_message( 'initializing by user', 'start' )
909	!
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	!
915	!-- Initialization will completely be done by the user
916	CALL user_init_3d_model
917
918	IF ( debug_output ) CALL debug_message( 'initializing by user', 'end' )
919
920	ENDIF
921
922	IF ( debug_output ) CALL debug_message( 'initializing statistics, boundary conditions, etc.', 'start' )
923
924	!
925	!-- Bottom boundary
926	IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2 ) THEN
927	u(nzb,:,:) = 0.0_wp
928	v(nzb,:,:) = 0.0_wp
929	ENDIF
930
931	!
932	!-- Apply channel flow boundary condition
933	IF ( TRIM( bc_uv_t ) == 'dirichlet_0' ) THEN
934	u(nzt+1,:,:) = 0.0_wp
935	v(nzt+1,:,:) = 0.0_wp
936	ENDIF
937
938	!
939	!-- Calculate virtual potential temperature
940	IF ( humidity ) vpt = pt * ( 1.0_wp + 0.61_wp * q )
941
942	!
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.
947	hom(:,1,5,:) = SPREAD( u(:,nys,nxl), 2, statistic_regions+1 )
948	hom(:,1,6,:) = SPREAD( v(:,nys,nxl), 2, statistic_regions+1 )
949	IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2) THEN
950	hom(nzb,1,5,:) = 0.0_wp
951	hom(nzb,1,6,:) = 0.0_wp
952	ENDIF
953	hom(:,1,7,:) = SPREAD( pt(:,nys,nxl), 2, statistic_regions+1 )
954
955	IF ( humidity ) THEN
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 )
961	!
962	!-- Store initial profile of mixing ratio and potential
963	!-- temperature
964	IF ( bulk_cloud_model .OR. cloud_droplets ) THEN
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
970	!
971	!-- Store initial scalar profile
972	IF ( passive_scalar ) THEN
973	hom(:,1,121,:) = SPREAD( s(:,nys,nxl), 2, statistic_regions+1 )
974	ENDIF
975
976	!
977	!-- Initialize the random number generators (from numerical recipes)
978	CALL random_function_ini
979
980	IF ( random_generator == 'random-parallel' ) THEN
981	CALL init_parallel_random_generator( nx, nys, nyn, nxl, nxr )
982	ENDIF
983	!
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
987	IF ( .NOT. humidity ) THEN
988	ref_state(:) = pt_reference
989	ELSE
990	ref_state(:) = vpt_reference
991	ENDIF
992	ELSE
993	IF ( .NOT. humidity ) THEN
994	ref_state(:) = pt_init(:)
995	ELSE
996	ref_state(:) = vpt(:,nys,nxl)
997	ENDIF
998	ENDIF
999
1000	!
1001	!-- For the moment, vertical velocity is zero
1002	w = 0.0_wp
1003
1004	!
1005	!-- Initialize array sums (must be defined in first call of pres)
1006	sums = 0.0_wp
1007
1008	!
1009	!-- In case of iterative solvers, p must get an initial value
1010	IF ( psolver(1:9) == 'multigrid' .OR. psolver == 'sor' ) p = 0.0_wp
1011	!
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	!
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
1022	CALL init_pt_anomaly
1023	ENDIF
1024
1025	!
1026	!-- If required, change the surface temperature at the start of the 3D run
1027	IF ( pt_surface_initial_change /= 0.0_wp ) THEN
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
1034	IF ( humidity .AND. q_surface_initial_change /= 0.0_wp ) &
1035	q(nzb,:,:) = q(nzb,:,:) + q_surface_initial_change
1036
1037	IF ( passive_scalar .AND. s_surface_initial_change /= 0.0_wp ) &
1038	s(nzb,:,:) = s(nzb,:,:) + s_surface_initial_change
1039
1040
1041	!
1042	!-- Initialize old and new time levels.
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
1045
1046	IF ( humidity ) THEN
1047	tq_m = 0.0_wp
1048	q_p = q
1049	ENDIF
1050
1051	IF ( passive_scalar ) THEN
1052	ts_m = 0.0_wp
1053	s_p = s
1054	ENDIF
1055
1056	IF ( debug_output ) CALL debug_message( 'initializing statistics, boundary conditions, etc.', 'end' )
1057
1058	ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data' .OR. &
1059	TRIM( initializing_actions ) == 'cyclic_fill' ) &
1060	THEN
1061
1062	IF ( debug_output ) CALL debug_message( 'initializing in case of restart / cyclic_fill', 'start' )
1063	!
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	!
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
1078	IF ( TRIM( initializing_actions ) == 'cyclic_fill' ) THEN
1079
1080	DO i = 0, io_blocks-1
1081	IF ( i == io_group ) THEN
1082	CALL rrd_read_parts_of_global
1083	ENDIF
1084	#if defined( __parallel )
1085	CALL MPI_BARRIER( comm2d, ierr )
1086	#endif
1087	ENDDO
1088
1089	ENDIF
1090
1091	!
1092	!-- Read processor specific binary data from restart file
1093	DO i = 0, io_blocks-1
1094	IF ( i == io_group ) THEN
1095	CALL rrd_local
1096	ENDIF
1097	#if defined( __parallel )
1098	CALL MPI_BARRIER( comm2d, ierr )
1099	#endif
1100	ENDDO
1101	!
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
1106
1107	!
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
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)
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	!
1132	!-- Initialization of the turbulence recycling method
1133	IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND. &
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
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.
1144	ALLOCATE( mean_inflow_profiles(nzb:nzt+1,1:num_mean_inflow_profiles) )
1145
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
1150	mean_inflow_profiles(:,1) = hom_sum(:,1,0) ! u
1151	u_init(:) = hom_sum(:,1,0)
1152	mean_inflow_profiles(:,2) = hom_sum(:,2,0) ! v
1153	v_init(:) = hom_sum(:,2,0)
1154	ENDIF
1155	mean_inflow_profiles(:,4) = hom_sum(:,4,0) ! pt
1156	pt_init(:) = hom_sum(:,4,0)
1157	IF ( humidity ) &
1158	mean_inflow_profiles(:,6) = hom_sum(:,41,0) ! q
1159	q_init(:) = hom_sum(:,41,0)
1160	IF ( passive_scalar ) &
1161	mean_inflow_profiles(:,7) = hom_sum(:,115,0) ! s
1162	s_init(:) = hom_sum(:,115,0)
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
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)
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
1178
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
1205	!
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
1210	DO j = nysg, nyng
1211	DO k = nzb, nzt+1
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)
1214	ENDDO
1215	ENDDO
1216	ENDDO
1217	ENDIF
1218
1219	!
1220	!-- Use these mean profiles at the inflow (provided that Dirichlet
1221	!-- conditions are used)
1222	IF ( bc_dirichlet_l ) THEN
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)
1227	w(k,j,nxlg:-1) = 0.0_wp
1228	pt(k,j,nxlg:-1) = mean_inflow_profiles(k,4)
1229	IF ( humidity ) &
1230	q(k,j,nxlg:-1) = mean_inflow_profiles(k,6)
1231	IF ( passive_scalar ) &
1232	s(k,j,nxlg:-1) = mean_inflow_profiles(k,7)
1233	ENDDO
1234	ENDDO
1235	ENDIF
1236
1237	!
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.
1242	IF ( inflow_damping_height == 9999999.9_wp ) THEN
1243	!
1244	!-- Default: use the inversion height calculated by the prerun; if
1245	!-- this is zero, inflow_damping_height must be explicitly
1246	!-- specified.
1247	IF ( hom_sum(nzb+6,pr_palm,0) /= 0.0_wp ) THEN
1248	inflow_damping_height = hom_sum(nzb+6,pr_palm,0)
1249	ELSE
1250	WRITE( message_string, * ) 'inflow_damping_height must be ', &
1251	'explicitly specified because&the inversion height ', &
1252	'calculated by the prerun is zero.'
1253	CALL message( 'init_3d_model', 'PA0318', 1, 2, 0, 6, 0 )
1254	ENDIF
1255
1256	ENDIF
1257
1258	IF ( inflow_damping_width == 9999999.9_wp ) THEN
1259	!
1260	!-- Default for the transition range: one tenth of the undamped
1261	!-- layer
1262	inflow_damping_width = 0.1_wp * inflow_damping_height
1263
1264	ENDIF
1265
1266	ALLOCATE( inflow_damping_factor(nzb:nzt+1) )
1267
1268	DO k = nzb, nzt+1
1269
1270	IF ( zu(k) <= inflow_damping_height ) THEN
1271	inflow_damping_factor(k) = 1.0_wp
1272	ELSEIF ( zu(k) <= ( inflow_damping_height + inflow_damping_width ) ) THEN
1273	inflow_damping_factor(k) = 1.0_wp - &
1274	( zu(k) - inflow_damping_height ) / &
1275	inflow_damping_width
1276	ELSE
1277	inflow_damping_factor(k) = 0.0_wp
1278	ENDIF
1279
1280	ENDDO
1281
1282	ENDIF
1283
1284	!
1285	!-- Inside buildings set velocities back to zero
1286	IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND. &
1287	topography /= 'flat' ) THEN
1288	!
1289	!-- Inside buildings set velocities back to zero.
1290	!-- Other scalars (pt, q, s, p, sa, ...) are ignored at present,
1291	!-- maybe revise later.
1292	DO i = nxlg, nxrg
1293	DO j = nysg, nyng
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
1302	ENDDO
1303	ENDDO
1304
1305	ENDIF
1306
1307	!
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)
1315	pt_p = pt; u_p = u; v_p = v; w_p = w
1316	IF ( humidity ) THEN
1317	q_p = q
1318	ENDIF
1319	IF ( passive_scalar ) s_p = s
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.
1324	tpt_m = 0.0_wp; tu_m = 0.0_wp; tv_m = 0.0_wp; tw_m = 0.0_wp
1325	IF ( humidity ) THEN
1326	tq_m = 0.0_wp
1327	ENDIF
1328	IF ( passive_scalar ) ts_m = 0.0_wp
1329
1330	IF ( debug_output ) CALL debug_message( 'initializing in case of restart / cyclic_fill', 'end' )
1331
1332	ELSE
1333	!
1334	!-- Actually this part of the programm should not be reached
1335	message_string = 'unknown initializing problem'
1336	CALL message( 'init_3d_model', 'PA0193', 1, 2, 0, 6, 0 )
1337	ENDIF
1338
1339
1340	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
1341	!
1342	!-- Initialize old timelevels needed for radiation boundary conditions
1343	IF ( bc_radiation_l ) THEN
1344	u_m_l(:,:,:) = u(:,:,1:2)
1345	v_m_l(:,:,:) = v(:,:,0:1)
1346	w_m_l(:,:,:) = w(:,:,0:1)
1347	ENDIF
1348	IF ( bc_radiation_r ) THEN
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
1353	IF ( bc_radiation_s ) THEN
1354	u_m_s(:,:,:) = u(:,0:1,:)
1355	v_m_s(:,:,:) = v(:,1:2,:)
1356	w_m_s(:,:,:) = w(:,0:1,:)
1357	ENDIF
1358	IF ( bc_radiation_n ) THEN
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
1363
1364	ENDIF
1365
1366	!
1367	!-- Calculate the initial volume flow at the right and north boundary
1368	IF ( conserve_volume_flow ) THEN
1369
1370	IF ( use_prescribed_profile_data ) THEN
1371
1372	volume_flow_initial_l = 0.0_wp
1373	volume_flow_area_l = 0.0_wp
1374
1375	IF ( nxr == nx ) THEN
1376	DO j = nys, nyn
1377	DO k = nzb+1, nzt
1378	volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
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	)
1388	ENDDO
1389	ENDDO
1390	ENDIF
1391
1392	IF ( nyn == ny ) THEN
1393	DO i = nxl, nxr
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	)
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
1421	volume_flow_initial_l = 0.0_wp
1422	volume_flow_area_l = 0.0_wp
1423
1424	IF ( nxr == nx ) THEN
1425	DO j = nys, nyn
1426	DO k = nzb+1, nzt
1427	volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
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	)
1436	ENDDO
1437	ENDDO
1438	ENDIF
1439
1440	IF ( nyn == ny ) THEN
1441	DO i = nxl, nxr
1442	DO k = nzb+1, nzt
1443	volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
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	)
1452	ENDDO
1453	ENDDO
1454	ENDIF
1455
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
1467	ELSEIF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
1468
1469	volume_flow_initial_l = 0.0_wp
1470	volume_flow_area_l = 0.0_wp
1471
1472	IF ( nxr == nx ) THEN
1473	DO j = nys, nyn
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	)
1484	ENDDO
1485	ENDDO
1486	ENDIF
1487
1488	IF ( nyn == ny ) THEN
1489	DO i = nxl, nxr
1490	DO k = nzb+1, nzt
1491	volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
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	)
1500	ENDDO
1501	ENDDO
1502	ENDIF
1503
1504	#if defined( __parallel )
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 )
1509
1510	#else
1511	volume_flow_initial = volume_flow_initial_l
1512	volume_flow_area = volume_flow_area_l
1513	#endif
1514
1515	ENDIF
1516
1517	!
1518	!-- In case of 'bulk_velocity' mode, volume_flow_initial is calculated
1519	!-- from u\|v_bulk instead
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
1524
1525	ENDIF
1526	!
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	!
1535	!-- Allocate memory for possible static input
1536	ALLOCATE( tmp_2d%var(nys:nyn,nxl:nxr) )
1537	tmp_2d%var = 0.0_wp
1538	!
1539	!-- Open the static input file
1540	#if defined( __netcdf )
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 )
1595	#endif
1596	DEALLOCATE( tmp_2d%var )
1597	ENDIF
1598	!
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.
1605	IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. &
1606	TRIM( initializing_actions ) /= 'cyclic_fill' ) THEN
1607
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
1615
1616	!
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.
1620	!-- ngp_2dh: number of grid points of a horizontal cross section through the
1621	!-- respective statistic region
1622	!-- ngp_3d: number of grid points of the respective statistic region
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
1693
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	!
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	!
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	!
1749	!-- Initialize quantities for special advections schemes
1750	CALL init_advec
1751
1752	!
1753	!-- Impose random perturbation on the horizontal velocity field and then
1754	!-- remove the divergences from the velocity field at the initial stage
1755	IF ( create_disturbances .AND. disturbance_energy_limit /= 0.0_wp .AND. &
1756	TRIM( initializing_actions ) /= 'read_restart_data' .AND. &
1757	TRIM( initializing_actions ) /= 'cyclic_fill' ) THEN
1758
1759	IF ( debug_output ) CALL debug_message( 'creating disturbances + applying pressure solver', 'start' )
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
1767	CALL disturb_field( 'u', tend, u )
1768	CALL disturb_field( 'v', tend, v )
1769
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
1785	n_sor = nsor_ini
1786	CALL pres
1787	n_sor = nsor
1788
1789	!$ACC END DATA
1790	!$ACC END DATA
1791
1792	IF ( debug_output ) CALL debug_message( 'creating disturbances + applying pressure solver', 'end' )
1793
1794	ENDIF
1795
1796	IF ( .NOT. ocean_mode ) THEN
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) )
1802	!
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)
1808	CALL message( 'init_3d_model', 'PA0142', 1, 2, 0, 6, 0 )
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
1822	ENDIF
1823
1824	!
1825	!-- If required, initialize particles
1826	IF ( agents_active ) CALL mas_init
1827	!
1828	!-- Initialization of synthetic turbulence generator
1829	IF ( use_syn_turb_gen ) CALL stg_init
1830	!
1831	!-- Initializing actions for all other modules
1832	CALL module_interface_init
1833	!
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
1837	!
1838	!-- Initialize surface data output
1839	IF ( surface_output ) CALL surface_data_output_init
1840	!
1841	!-- Initialize the ws-scheme.
1842	IF ( ws_scheme_sca .OR. ws_scheme_mom ) CALL ws_init
1843	!
1844	!-- Perform post-initializing checks for all other modules
1845	CALL module_interface_init_checks
1846
1847	!
1848	!-- Setting weighting factors for calculation of perturbation pressure
1849	!-- and turbulent quantities from the RK substeps
1850	IF ( TRIM(timestep_scheme) == 'runge-kutta-3' ) THEN ! for RK3-method
1851
1852	weight_substep(1) = 1._wp/6._wp
1853	weight_substep(2) = 3._wp/10._wp
1854	weight_substep(3) = 8._wp/15._wp
1855
1856	weight_pres(1) = 1._wp/3._wp
1857	weight_pres(2) = 5._wp/12._wp
1858	weight_pres(3) = 1._wp/4._wp
1859
1860	ELSEIF ( TRIM(timestep_scheme) == 'runge-kutta-2' ) THEN ! for RK2-method
1861
1862	weight_substep(1) = 1._wp/2._wp
1863	weight_substep(2) = 1._wp/2._wp
1864
1865	weight_pres(1) = 1._wp/2._wp
1866	weight_pres(2) = 1._wp/2._wp
1867
1868	ELSE ! for Euler-method
1869
1870	weight_substep(1) = 1.0_wp
1871	weight_pres(1) = 1.0_wp
1872
1873	ENDIF
1874
1875	!
1876	!-- Initialize Rayleigh damping factors
1877	rdf = 0.0_wp
1878	rdf_sc = 0.0_wp
1879	IF ( rayleigh_damping_factor /= 0.0_wp ) THEN
1880
1881	IF ( .NOT. ocean_mode ) THEN
1882	DO k = nzb+1, nzt
1883	IF ( zu(k) >= rayleigh_damping_height ) THEN
1884	rdf(k) = rayleigh_damping_factor * &
1885	( SIN( pi * 0.5_wp * ( zu(k) - rayleigh_damping_height ) &
1886	/ ( zu(nzt) - rayleigh_damping_height ) ) &
1887	)**2
1888	ENDIF
1889	ENDDO
1890	ELSE
1891	!
1892	!-- In ocean mode, rayleigh damping is applied in the lower part of the
1893	!-- model domain
1894	DO k = nzt, nzb+1, -1
1895	IF ( zu(k) <= rayleigh_damping_height ) THEN
1896	rdf(k) = rayleigh_damping_factor * &
1897	( SIN( pi * 0.5_wp * ( rayleigh_damping_height - zu(k) ) &
1898	/ ( rayleigh_damping_height - zu(nzb+1) ) ) &
1899	)**2
1900	ENDIF
1901	ENDDO
1902	ENDIF
1903
1904	ENDIF
1905	IF ( scalar_rayleigh_damping ) rdf_sc = rdf
1906
1907	!
1908	!-- Initialize the starting level and the vertical smoothing factor used for
1909	!-- the external pressure gradient
1910	dp_smooth_factor = 1.0_wp
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
1921	dp_smooth_factor(:dp_level_ind_b) = 0.0_wp
1922	DO k = dp_level_ind_b+1, nzt
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 ) ) )
1926	ENDDO
1927	ENDIF
1928	ENDIF
1929
1930	!
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.
1934	ptdf_x = 0.0_wp
1935	ptdf_y = 0.0_wp
1936	IF ( bc_lr_dirrad ) THEN
1937	DO i = nxl, nxr
1938	IF ( ( i * dx ) < pt_damping_width ) THEN
1939	ptdf_x(i) = pt_damping_factor * ( SIN( pi * 0.5_wp * &
1940	REAL( pt_damping_width - i * dx, KIND=wp ) / ( &
1941	REAL( pt_damping_width, KIND=wp ) ) ) )**2
1942	ENDIF
1943	ENDDO
1944	ELSEIF ( bc_lr_raddir ) THEN
1945	DO i = nxl, nxr
1946	IF ( ( i * dx ) > ( nx * dx - pt_damping_width ) ) THEN
1947	ptdf_x(i) = pt_damping_factor * &
1948	SIN( pi * 0.5_wp * &
1949	( ( i - nx ) * dx + pt_damping_width ) / &
1950	REAL( pt_damping_width, KIND=wp ) )**2
1951	ENDIF
1952	ENDDO
1953	ELSEIF ( bc_ns_dirrad ) THEN
1954	DO j = nys, nyn
1955	IF ( ( j * dy ) > ( ny * dy - pt_damping_width ) ) THEN
1956	ptdf_y(j) = pt_damping_factor * &
1957	SIN( pi * 0.5_wp * &
1958	( ( j - ny ) * dy + pt_damping_width ) / &
1959	REAL( pt_damping_width, KIND=wp ) )**2
1960	ENDIF
1961	ENDDO
1962	ELSEIF ( bc_ns_raddir ) THEN
1963	DO j = nys, nyn
1964	IF ( ( j * dy ) < pt_damping_width ) THEN
1965	ptdf_y(j) = pt_damping_factor * &
1966	SIN( pi * 0.5_wp * &
1967	( pt_damping_width - j * dy ) / &
1968	REAL( pt_damping_width, KIND=wp ) )**2
1969	ENDIF
1970	ENDDO
1971	ENDIF
1972
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 )
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
1983
1984
1985	CALL location_message( 'model initialization', 'finished' )
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: