Home

Context Navigation

source: palm/trunk/SOURCE/nesting_offl_mod.f90 @ 4230

Last change on this file since 4230 was 4230, checked in by suehring, 5 years ago
Several bugfixes: profile initialization of chemical species in restart runs; Runge-Kutta tendency array not initialized in chemistry model in restart runs; fixed determination of time indices for chemical emissions (introduced with commit -4227); Update chemistry profiles in offline nesting; initialize canopy resistances for greened building walls (even if green fraction is zero)
Property svn:keywords set to `Id`
File size: 110.5 KB

Line
1	!> @file nesting_offl_mod.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: nesting_offl_mod.f90 4230 2019-09-11 13:58:14Z suehring $
27	! Update mean chemistry profiles. These are also used for rayleigh damping.
28	!
29	! 4227 2019-09-10 18:04:34Z gronemeier
30	! implement new palm_date_time_mod
31	!
32	! - Data input moved into nesting_offl_mod
33	! - check rephrased
34	! - time variable is now relative to time_utc_init
35	! - Define module specific data type for offline nesting in nesting_offl_mod
36	!
37	! 4182 2019-08-22 15:20:23Z scharf
38	! Corrected "Former revisions" section
39	!
40	! 4169 2019-08-19 13:54:35Z suehring
41	! Additional check added.
42	!
43	! 4168 2019-08-16 13:50:17Z suehring
44	! Replace function get_topography_top_index by topo_top_ind
45	!
46	! 4125 2019-07-29 13:31:44Z suehring
47	! In order to enable netcdf parallel access, allocate dummy arrays for the
48	! lateral boundary data on cores that actually do not belong to these
49	! boundaries.
50	!
51	! 4079 2019-07-09 18:04:41Z suehring
52	! - Set boundary condition for w at nzt+1 at the lateral boundaries, even
53	! though these won't enter the numerical solution. However, due to the mass
54	! conservation these values might some up to very large values which will
55	! occur in the run-control file
56	! - Bugfix in offline nesting of chemical species
57	! - Do not set Neumann conditions for TKE and passive scalar
58	!
59	! 4022 2019-06-12 11:52:39Z suehring
60	! Detection of boundary-layer depth in stable boundary layer on basis of
61	! boundary data improved
62	! Routine for boundary-layer depth calculation renamed and made public
63	!
64	! 3987 2019-05-22 09:52:13Z kanani
65	! Introduce alternative switch for debug output during timestepping
66	!
67	! 3964 2019-05-09 09:48:32Z suehring
68	! Ensure that veloctiy term in calculation of bulk Richardson number does not
69	! become zero
70	!
71	! 3937 2019-04-29 15:09:07Z suehring
72	! Set boundary conditon on upper-left and upper-south grid point for the u- and
73	! v-component, respectively.
74	!
75	! 3891 2019-04-12 17:52:01Z suehring
76	! Bugfix, do not overwrite lateral and top boundary data in case of restart
77	! runs.
78	!
79	! 3885 2019-04-11 11:29:34Z kanani
80	! Changes related to global restructuring of location messages and introduction
81	! of additional debug messages
82	!
83	!
84	! Do local data exchange for chemistry variables only when boundary data is
85	! coming from dynamic file
86	!
87	! 3737 2019-02-12 16:57:06Z suehring
88	! Introduce mesoscale nesting for chemical species
89	!
90	! 3705 2019-01-29 19:56:39Z suehring
91	! Formatting adjustments
92	!
93	! 3704 2019-01-29 19:51:41Z suehring
94	! Check implemented for offline nesting in child domain
95	!
96	! Initial Revision:
97	! - separate offline nesting from large_scale_nudging_mod
98	! - revise offline nesting, adjust for usage of synthetic turbulence generator
99	! - adjust Rayleigh damping depending on the time-depending atmospheric
100	! conditions
101	!
102	!
103	! Description:
104	! ------------
105	!> Offline nesting in larger-scale models. Boundary conditions for the simulation
106	!> are read from NetCDF file and are prescribed onto the respective arrays.
107	!> Further, a mass-flux correction is performed to maintain the mass balance.
108	!--------------------------------------------------------------------------------!
109	MODULE nesting_offl_mod
110
111	USE arrays_3d, &
112	ONLY: dzw, &
113	e, &
114	diss, &
115	pt, &
116	pt_init, &
117	q, &
118	q_init, &
119	rdf, &
120	rdf_sc, &
121	s, &
122	u, &
123	u_init, &
124	ug, &
125	v, &
126	v_init, &
127	vg, &
128	w, &
129	zu, &
130	zw
131
132	USE basic_constants_and_equations_mod, &
133	ONLY: g, &
134	pi
135
136	USE chem_modules, &
137	ONLY: chem_species
138
139	USE control_parameters, &
140	ONLY: air_chemistry, &
141	bc_dirichlet_l, &
142	bc_dirichlet_n, &
143	bc_dirichlet_r, &
144	bc_dirichlet_s, &
145	coupling_char, &
146	dt_3d, &
147	dz, &
148	constant_diffusion, &
149	child_domain, &
150	debug_output_timestep, &
151	end_time, &
152	humidity, &
153	initializing_actions, &
154	message_string, &
155	nesting_offline, &
156	neutral, &
157	passive_scalar, &
158	rans_mode, &
159	rans_tke_e, &
160	rayleigh_damping_factor, &
161	rayleigh_damping_height, &
162	spinup_time, &
163	time_since_reference_point, &
164	volume_flow
165
166	USE cpulog, &
167	ONLY: cpu_log, &
168	log_point, &
169	log_point_s
170
171	USE grid_variables
172
173	USE indices, &
174	ONLY: nbgp, nx, nxl, nxlg, nxlu, nxr, nxrg, ny, nys, &
175	nysv, nysg, nyn, nyng, nzb, nz, nzt, &
176	topo_top_ind, &
177	wall_flags_0
178
179	USE kinds
180
181	USE netcdf_data_input_mod, &
182	ONLY: check_existence, &
183	close_input_file, &
184	get_dimension_length, &
185	get_variable, &
186	get_variable_pr, &
187	input_pids_dynamic, &
188	inquire_num_variables, &
189	inquire_variable_names, &
190	input_file_dynamic, &
191	num_var_pids, &
192	open_read_file, &
193	pids_id
194
195	USE palm_date_time_mod, &
196	ONLY: get_date_time
197
198	USE pegrid
199
200	IMPLICIT NONE
201
202	!
203	!-- Define data type for nesting in larger-scale models like COSMO.
204	!-- Data type comprises u, v, w, pt, and q at lateral and top boundaries.
205	TYPE nest_offl_type
206
207	CHARACTER(LEN=16) :: char_l = 'ls_forcing_left_' !< leading substring for variables at left boundary
208	CHARACTER(LEN=17) :: char_n = 'ls_forcing_north_' !< leading substring for variables at north boundary
209	CHARACTER(LEN=17) :: char_r = 'ls_forcing_right_' !< leading substring for variables at right boundary
210	CHARACTER(LEN=17) :: char_s = 'ls_forcing_south_' !< leading substring for variables at south boundary
211	CHARACTER(LEN=15) :: char_t = 'ls_forcing_top_' !< leading substring for variables at top boundary
212
213	CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE :: var_names !< list of variable in dynamic input file
214	CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE :: var_names_chem_l !< names of mesoscale nested chemistry variables at left boundary
215	CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE :: var_names_chem_n !< names of mesoscale nested chemistry variables at north boundary
216	CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE :: var_names_chem_r !< names of mesoscale nested chemistry variables at right boundary
217	CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE :: var_names_chem_s !< names of mesoscale nested chemistry variables at south boundary
218	CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE :: var_names_chem_t !< names of mesoscale nested chemistry variables at top boundary
219
220	INTEGER(iwp) :: nt !< number of time levels in dynamic input file
221	INTEGER(iwp) :: nzu !< number of vertical levels on scalar grid in dynamic input file
222	INTEGER(iwp) :: nzw !< number of vertical levels on w grid in dynamic input file
223	INTEGER(iwp) :: tind !< time index for reference time in mesoscale-offline nesting
224	INTEGER(iwp) :: tind_p !< time index for following time in mesoscale-offline nesting
225
226	LOGICAL :: init = .FALSE. !< flag indicating that offline nesting is already initialized
227
228	LOGICAL, DIMENSION(:), ALLOCATABLE :: chem_from_file_l !< flags inidicating whether left boundary data for chemistry is in dynamic input file
229	LOGICAL, DIMENSION(:), ALLOCATABLE :: chem_from_file_n !< flags inidicating whether north boundary data for chemistry is in dynamic input file
230	LOGICAL, DIMENSION(:), ALLOCATABLE :: chem_from_file_r !< flags inidicating whether right boundary data for chemistry is in dynamic input file
231	LOGICAL, DIMENSION(:), ALLOCATABLE :: chem_from_file_s !< flags inidicating whether south boundary data for chemistry is in dynamic input file
232	LOGICAL, DIMENSION(:), ALLOCATABLE :: chem_from_file_t !< flags inidicating whether top boundary data for chemistry is in dynamic input file
233
234	REAL(wp), DIMENSION(:), ALLOCATABLE :: surface_pressure !< time dependent surface pressure
235	REAL(wp), DIMENSION(:), ALLOCATABLE :: time !< time levels in dynamic input file
236	REAL(wp), DIMENSION(:), ALLOCATABLE :: zu_atmos !< vertical levels at scalar grid in dynamic input file
237	REAL(wp), DIMENSION(:), ALLOCATABLE :: zw_atmos !< vertical levels at w grid in dynamic input file
238
239	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ug !< domain-averaged geostrophic component
240	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: vg !< domain-averaged geostrophic component
241
242	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: u_left !< u-component at left boundary
243	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: v_left !< v-component at left boundary
244	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: w_left !< w-component at left boundary
245	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: q_left !< mixing ratio at left boundary
246	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pt_left !< potentital temperautre at left boundary
247
248	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: u_north !< u-component at north boundary
249	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: v_north !< v-component at north boundary
250	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: w_north !< w-component at north boundary
251	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: q_north !< mixing ratio at north boundary
252	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pt_north !< potentital temperautre at north boundary
253
254	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: u_right !< u-component at right boundary
255	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: v_right !< v-component at right boundary
256	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: w_right !< w-component at right boundary
257	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: q_right !< mixing ratio at right boundary
258	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pt_right !< potentital temperautre at right boundary
259
260	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: u_south !< u-component at south boundary
261	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: v_south !< v-component at south boundary
262	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: w_south !< w-component at south boundary
263	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: q_south !< mixing ratio at south boundary
264	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pt_south !< potentital temperautre at south boundary
265
266	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: u_top !< u-component at top boundary
267	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: v_top !< v-component at top boundary
268	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: w_top !< w-component at top boundary
269	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: q_top !< mixing ratio at top boundary
270	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pt_top !< potentital temperautre at top boundary
271
272	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: chem_left !< chemical species at left boundary
273	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: chem_north !< chemical species at left boundary
274	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: chem_right !< chemical species at left boundary
275	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: chem_south !< chemical species at left boundary
276	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: chem_top !< chemical species at left boundary
277
278	END TYPE nest_offl_type
279
280	REAL(wp) :: fac_dt !< interpolation factor
281	REAL(wp) :: time_utc_init !< time in seconds-of-day of origin_date_time
282	REAL(wp) :: zi_ribulk = 0.0_wp !< boundary-layer depth according to bulk Richardson criterion, i.e. the height where Ri_bulk exceeds the critical
283	!< bulk Richardson number of 0.2
284
285	TYPE(nest_offl_type) :: nest_offl !< data structure for data input at lateral and top boundaries (provided by Inifor)
286
287	SAVE
288	PRIVATE
289	!
290	!-- Public subroutines
291	PUBLIC nesting_offl_bc, &
292	nesting_offl_calc_zi, &
293	nesting_offl_check_parameters, &
294	nesting_offl_geostrophic_wind, &
295	nesting_offl_header, &
296	nesting_offl_init, &
297	nesting_offl_input, &
298	nesting_offl_interpolation_factor, &
299	nesting_offl_mass_conservation, &
300	nesting_offl_parin
301	!
302	!-- Public variables
303	PUBLIC zi_ribulk
304
305	INTERFACE nesting_offl_bc
306	MODULE PROCEDURE nesting_offl_bc
307	END INTERFACE nesting_offl_bc
308
309	INTERFACE nesting_offl_calc_zi
310	MODULE PROCEDURE nesting_offl_calc_zi
311	END INTERFACE nesting_offl_calc_zi
312
313	INTERFACE nesting_offl_check_parameters
314	MODULE PROCEDURE nesting_offl_check_parameters
315	END INTERFACE nesting_offl_check_parameters
316
317	INTERFACE nesting_offl_geostrophic_wind
318	MODULE PROCEDURE nesting_offl_geostrophic_wind
319	END INTERFACE nesting_offl_geostrophic_wind
320
321	INTERFACE nesting_offl_header
322	MODULE PROCEDURE nesting_offl_header
323	END INTERFACE nesting_offl_header
324
325	INTERFACE nesting_offl_init
326	MODULE PROCEDURE nesting_offl_init
327	END INTERFACE nesting_offl_init
328
329	INTERFACE nesting_offl_input
330	MODULE PROCEDURE nesting_offl_input
331	END INTERFACE nesting_offl_input
332
333	INTERFACE nesting_offl_interpolation_factor
334	MODULE PROCEDURE nesting_offl_interpolation_factor
335	END INTERFACE nesting_offl_interpolation_factor
336
337	INTERFACE nesting_offl_mass_conservation
338	MODULE PROCEDURE nesting_offl_mass_conservation
339	END INTERFACE nesting_offl_mass_conservation
340
341	INTERFACE nesting_offl_parin
342	MODULE PROCEDURE nesting_offl_parin
343	END INTERFACE nesting_offl_parin
344
345	CONTAINS
346
347	!------------------------------------------------------------------------------!
348	! Description:
349	! ------------
350	!> Reads data at lateral and top boundaries derived from larger-scale model.
351	!------------------------------------------------------------------------------!
352	SUBROUTINE nesting_offl_input
353
354	INTEGER(iwp) :: n !< running index for chemistry variables
355	INTEGER(iwp) :: t !< running index time dimension
356
357	!
358	!-- Initialize INIFOR forcing in first call.
359	IF ( .NOT. nest_offl%init ) THEN
360	#if defined ( __netcdf )
361	!
362	!-- Open file in read-only mode
363	CALL open_read_file( TRIM( input_file_dynamic ) // &
364	TRIM( coupling_char ), pids_id )
365	!
366	!-- At first, inquire all variable names.
367	CALL inquire_num_variables( pids_id, num_var_pids )
368	!
369	!-- Allocate memory to store variable names.
370	ALLOCATE( nest_offl%var_names(1:num_var_pids) )
371	CALL inquire_variable_names( pids_id, nest_offl%var_names )
372	!
373	!-- Read time dimension, allocate memory and finally read time array
374	CALL get_dimension_length( pids_id, nest_offl%nt, 'time' )
375
376	IF ( check_existence( nest_offl%var_names, 'time' ) ) THEN
377	ALLOCATE( nest_offl%time(0:nest_offl%nt-1) )
378	CALL get_variable( pids_id, 'time', nest_offl%time )
379	ENDIF
380	!
381	!-- Read vertical dimension of scalar und w grid
382	CALL get_dimension_length( pids_id, nest_offl%nzu, 'z' )
383	CALL get_dimension_length( pids_id, nest_offl%nzw, 'zw' )
384
385	IF ( check_existence( nest_offl%var_names, 'z' ) ) THEN
386	ALLOCATE( nest_offl%zu_atmos(1:nest_offl%nzu) )
387	CALL get_variable( pids_id, 'z', nest_offl%zu_atmos )
388	ENDIF
389	IF ( check_existence( nest_offl%var_names, 'zw' ) ) THEN
390	ALLOCATE( nest_offl%zw_atmos(1:nest_offl%nzw) )
391	CALL get_variable( pids_id, 'zw', nest_offl%zw_atmos )
392	ENDIF
393	!
394	!-- Read surface pressure
395	IF ( check_existence( nest_offl%var_names, &
396	'surface_forcing_surface_pressure' ) ) THEN
397	ALLOCATE( nest_offl%surface_pressure(0:nest_offl%nt-1) )
398	CALL get_variable( pids_id, &
399	'surface_forcing_surface_pressure', &
400	nest_offl%surface_pressure )
401	ENDIF
402	!
403	!-- Close input file
404	CALL close_input_file( pids_id )
405	#endif
406	ENDIF
407	!
408	!-- Check if dynamic driver data input is required.
409	IF ( nest_offl%time(nest_offl%tind_p) <= &
410	MAX( time_since_reference_point, 0.0_wp) + time_utc_init .OR. &
411	.NOT. nest_offl%init ) THEN
412	CONTINUE
413	!
414	!-- Return otherwise
415	ELSE
416	RETURN
417	ENDIF
418	!
419	!-- CPU measurement
420	CALL cpu_log( log_point_s(86), 'NetCDF input forcing', 'start' )
421
422	!
423	!-- Obtain time index for current point in time. Note, the time coordinate
424	!-- in the input file is relative to time_utc_init. Since time_since_...
425	!-- is negativ when spinup is used, use MAX function to obtain correct
426	!-- time at the beginning.
427	nest_offl%tind = MINLOC( ABS( nest_offl%time - ( &
428	time_utc_init + &
429	MAX( time_since_reference_point, 0.0_wp) )&
430	), DIM = 1 ) - 1
431	nest_offl%tind_p = nest_offl%tind + 1
432	!
433	!-- Open file in read-only mode
434	#if defined ( __netcdf )
435	CALL open_read_file( TRIM( input_file_dynamic ) // &
436	TRIM( coupling_char ), pids_id )
437	!
438	!-- Read geostrophic wind components
439	DO t = nest_offl%tind, nest_offl%tind_p
440	CALL get_variable_pr( pids_id, 'ls_forcing_ug', t+1, &
441	nest_offl%ug(t-nest_offl%tind,nzb+1:nzt) )
442	CALL get_variable_pr( pids_id, 'ls_forcing_vg', t+1, &
443	nest_offl%vg(t-nest_offl%tind,nzb+1:nzt) )
444	ENDDO
445	!
446	!-- Read data at lateral and top boundaries. Please note, at left and
447	!-- right domain boundary, yz-layers are read for u, v, w, pt and q.
448	!-- For the v-component, the data starts at nysv, while for the other
449	!-- quantities the data starts at nys. This is equivalent at the north
450	!-- and south domain boundary for the u-component.
451	!-- Note, lateral data is also accessed by parallel IO, which is the reason
452	!-- why different arguments are passed depending on the boundary control
453	!-- flags. Cores that do not belong to the respective boundary just make
454	!-- a dummy read with count = 0, just in order to participate the collective
455	!-- operation.
456	!-- Read data for western boundary
457	CALL get_variable( pids_id, 'ls_forcing_left_u', &
458	nest_offl%u_left, & ! array to be read
459	MERGE( nys+1, 1, bc_dirichlet_l), & ! start index y direction
460	MERGE( nzb+1, 1, bc_dirichlet_l), & ! start index z direction
461	MERGE( nest_offl%tind+1, 1, bc_dirichlet_l), & ! start index time dimension
462	MERGE( nyn-nys+1, 0, bc_dirichlet_l), & ! number of elements along y
463	MERGE( nest_offl%nzu, 0, bc_dirichlet_l), & ! number of elements alogn z
464	MERGE( 2, 0, bc_dirichlet_l), & ! number of time steps (2 or 0)
465	.TRUE. ) ! parallel IO when compiled accordingly
466
467	CALL get_variable( pids_id, 'ls_forcing_left_v', &
468	nest_offl%v_left, &
469	MERGE( nysv, 1, bc_dirichlet_l), &
470	MERGE( nzb+1, 1, bc_dirichlet_l), &
471	MERGE( nest_offl%tind+1, 1, bc_dirichlet_l), &
472	MERGE( nyn-nysv+1, 0, bc_dirichlet_l), &
473	MERGE( nest_offl%nzu, 0, bc_dirichlet_l), &
474	MERGE( 2, 0, bc_dirichlet_l), &
475	.TRUE. )
476
477	CALL get_variable( pids_id, 'ls_forcing_left_w', &
478	nest_offl%w_left, &
479	MERGE( nys+1, 1, bc_dirichlet_l), &
480	MERGE( nzb+1, 1, bc_dirichlet_l), &
481	MERGE( nest_offl%tind+1, 1, bc_dirichlet_l), &
482	MERGE( nyn-nys+1, 0, bc_dirichlet_l), &
483	MERGE( nest_offl%nzw, 0, bc_dirichlet_l), &
484	MERGE( 2, 0, bc_dirichlet_l), &
485	.TRUE. )
486
487	IF ( .NOT. neutral ) THEN
488	CALL get_variable( pids_id, 'ls_forcing_left_pt', &
489	nest_offl%pt_left, &
490	MERGE( nys+1, 1, bc_dirichlet_l), &
491	MERGE( nzb+1, 1, bc_dirichlet_l), &
492	MERGE( nest_offl%tind+1, 1, bc_dirichlet_l), &
493	MERGE( nyn-nys+1, 0, bc_dirichlet_l), &
494	MERGE( nest_offl%nzu, 0, bc_dirichlet_l), &
495	MERGE( 2, 0, bc_dirichlet_l), &
496	.TRUE. )
497	ENDIF
498
499	IF ( humidity ) THEN
500	CALL get_variable( pids_id, 'ls_forcing_left_qv', &
501	nest_offl%q_left, &
502	MERGE( nys+1, 1, bc_dirichlet_l), &
503	MERGE( nzb+1, 1, bc_dirichlet_l), &
504	MERGE( nest_offl%tind+1, 1, bc_dirichlet_l), &
505	MERGE( nyn-nys+1, 0, bc_dirichlet_l), &
506	MERGE( nest_offl%nzu, 0, bc_dirichlet_l), &
507	MERGE( 2, 0, bc_dirichlet_l), &
508	.TRUE. )
509	ENDIF
510
511	IF ( air_chemistry ) THEN
512	DO n = 1, UBOUND(nest_offl%var_names_chem_l, 1)
513	IF ( check_existence( nest_offl%var_names, &
514	nest_offl%var_names_chem_l(n) ) ) THEN
515	CALL get_variable( pids_id, &
516	TRIM( nest_offl%var_names_chem_l(n) ), &
517	nest_offl%chem_left(:,:,:,n), &
518	MERGE( nys+1, 1, bc_dirichlet_l), &
519	MERGE( nzb+1, 1, bc_dirichlet_l), &
520	MERGE( nest_offl%tind+1, 1, bc_dirichlet_l), &
521	MERGE( nyn-nys+1, 0, bc_dirichlet_l), &
522	MERGE( nest_offl%nzu, 0, bc_dirichlet_l), &
523	MERGE( 2, 0, bc_dirichlet_l), &
524	.TRUE. )
525	nest_offl%chem_from_file_l(n) = .TRUE.
526	ENDIF
527	ENDDO
528	ENDIF
529	!
530	!-- Read data for eastern boundary
531	CALL get_variable( pids_id, 'ls_forcing_right_u', &
532	nest_offl%u_right, &
533	MERGE( nys+1, 1, bc_dirichlet_r), &
534	MERGE( nzb+1, 1, bc_dirichlet_r), &
535	MERGE( nest_offl%tind+1, 1, bc_dirichlet_r), &
536	MERGE( nyn-nys+1, 0, bc_dirichlet_r), &
537	MERGE( nest_offl%nzu, 0, bc_dirichlet_r), &
538	MERGE( 2, 0, bc_dirichlet_r), &
539	.TRUE. )
540
541	CALL get_variable( pids_id, 'ls_forcing_right_v', &
542	nest_offl%v_right, &
543	MERGE( nysv, 1, bc_dirichlet_r), &
544	MERGE( nzb+1, 1, bc_dirichlet_r), &
545	MERGE( nest_offl%tind+1, 1, bc_dirichlet_r), &
546	MERGE( nyn-nysv+1, 0, bc_dirichlet_r), &
547	MERGE( nest_offl%nzu, 0, bc_dirichlet_r), &
548	MERGE( 2, 0, bc_dirichlet_r), &
549	.TRUE. )
550
551	CALL get_variable( pids_id, 'ls_forcing_right_w', &
552	nest_offl%w_right, &
553	MERGE( nys+1, 1, bc_dirichlet_r), &
554	MERGE( nzb+1, 1, bc_dirichlet_r), &
555	MERGE( nest_offl%tind+1, 1, bc_dirichlet_r), &
556	MERGE( nyn-nys+1, 0, bc_dirichlet_r), &
557	MERGE( nest_offl%nzw, 0, bc_dirichlet_r), &
558	MERGE( 2, 0, bc_dirichlet_r), &
559	.TRUE. )
560
561	IF ( .NOT. neutral ) THEN
562	CALL get_variable( pids_id, 'ls_forcing_right_pt', &
563	nest_offl%pt_right, &
564	MERGE( nys+1, 1, bc_dirichlet_r), &
565	MERGE( nzb+1, 1, bc_dirichlet_r), &
566	MERGE( nest_offl%tind+1, 1, bc_dirichlet_r), &
567	MERGE( nyn-nys+1, 0, bc_dirichlet_r), &
568	MERGE( nest_offl%nzu, 0, bc_dirichlet_r), &
569	MERGE( 2, 0, bc_dirichlet_r), &
570	.TRUE. )
571	ENDIF
572
573	IF ( humidity ) THEN
574	CALL get_variable( pids_id, 'ls_forcing_right_qv', &
575	nest_offl%q_right, &
576	MERGE( nys+1, 1, bc_dirichlet_r), &
577	MERGE( nzb+1, 1, bc_dirichlet_r), &
578	MERGE( nest_offl%tind+1, 1, bc_dirichlet_r), &
579	MERGE( nyn-nys+1, 0, bc_dirichlet_r), &
580	MERGE( nest_offl%nzu, 0, bc_dirichlet_r), &
581	MERGE( 2, 0, bc_dirichlet_r), &
582	.TRUE. )
583	ENDIF
584
585	IF ( air_chemistry ) THEN
586	DO n = 1, UBOUND(nest_offl%var_names_chem_r, 1)
587	IF ( check_existence( nest_offl%var_names, &
588	nest_offl%var_names_chem_r(n) ) ) THEN
589	CALL get_variable( pids_id, &
590	TRIM( nest_offl%var_names_chem_r(n) ), &
591	nest_offl%chem_right(:,:,:,n), &
592	MERGE( nys+1, 1, bc_dirichlet_r), &
593	MERGE( nzb+1, 1, bc_dirichlet_r), &
594	MERGE( nest_offl%tind+1, 1, bc_dirichlet_r), &
595	MERGE( nyn-nys+1, 0, bc_dirichlet_r), &
596	MERGE( nest_offl%nzu, 0, bc_dirichlet_r), &
597	MERGE( 2, 0, bc_dirichlet_r), &
598	.TRUE. )
599	nest_offl%chem_from_file_r(n) = .TRUE.
600	ENDIF
601	ENDDO
602	ENDIF
603	!
604	!-- Read data for northern boundary
605	CALL get_variable( pids_id, 'ls_forcing_north_u', & ! array to be read
606	nest_offl%u_north, & ! start index x direction
607	MERGE( nxlu, 1, bc_dirichlet_n ), & ! start index z direction
608	MERGE( nzb+1, 1, bc_dirichlet_n ), & ! start index time dimension
609	MERGE( nest_offl%tind+1, 1, bc_dirichlet_n ), & ! number of elements along x
610	MERGE( nxr-nxlu+1, 0, bc_dirichlet_n ), & ! number of elements alogn z
611	MERGE( nest_offl%nzu, 0, bc_dirichlet_n ), & ! number of time steps (2 or 0)
612	MERGE( 2, 0, bc_dirichlet_n ), & ! parallel IO when compiled accordingly
613	.TRUE. )
614
615	CALL get_variable( pids_id, 'ls_forcing_north_v', & ! array to be read
616	nest_offl%v_north, & ! start index x direction
617	MERGE( nxl+1, 1, bc_dirichlet_n ), & ! start index z direction
618	MERGE( nzb+1, 1, bc_dirichlet_n ), & ! start index time dimension
619	MERGE( nest_offl%tind+1, 1, bc_dirichlet_n ), & ! number of elements along x
620	MERGE( nxr-nxl+1, 0, bc_dirichlet_n ), & ! number of elements alogn z
621	MERGE( nest_offl%nzu, 0, bc_dirichlet_n ), & ! number of time steps (2 or 0)
622	MERGE( 2, 0, bc_dirichlet_n ), & ! parallel IO when compiled accordingly
623	.TRUE. )
624
625	CALL get_variable( pids_id, 'ls_forcing_north_w', & ! array to be read
626	nest_offl%w_north, & ! start index x direction
627	MERGE( nxl+1, 1, bc_dirichlet_n ), & ! start index z direction
628	MERGE( nzb+1, 1, bc_dirichlet_n ), & ! start index time dimension
629	MERGE( nest_offl%tind+1, 1, bc_dirichlet_n ), & ! number of elements along x
630	MERGE( nxr-nxl+1, 0, bc_dirichlet_n ), & ! number of elements alogn z
631	MERGE( nest_offl%nzw, 0, bc_dirichlet_n ), & ! number of time steps (2 or 0)
632	MERGE( 2, 0, bc_dirichlet_n ), & ! parallel IO when compiled accordingly
633	.TRUE. )
634
635	IF ( .NOT. neutral ) THEN
636	CALL get_variable( pids_id, 'ls_forcing_north_pt', & ! array to be read
637	nest_offl%pt_north, & ! start index x direction
638	MERGE( nxl+1, 1, bc_dirichlet_n ), & ! start index z direction
639	MERGE( nzb+1, 1, bc_dirichlet_n ), & ! start index time dimension
640	MERGE( nest_offl%tind+1, 1, bc_dirichlet_n ), & ! number of elements along x
641	MERGE( nxr-nxl+1, 0, bc_dirichlet_n ), & ! number of elements alogn z
642	MERGE( nest_offl%nzu, 0, bc_dirichlet_n ), & ! number of time steps (2 or 0)
643	MERGE( 2, 0, bc_dirichlet_n ), & ! parallel IO when compiled accordingly
644	.TRUE. )
645	ENDIF
646	IF ( humidity ) THEN
647	CALL get_variable( pids_id, 'ls_forcing_north_qv', & ! array to be read
648	nest_offl%q_north, & ! start index x direction
649	MERGE( nxl+1, 1, bc_dirichlet_n ), & ! start index z direction
650	MERGE( nzb+1, 1, bc_dirichlet_n ), & ! start index time dimension
651	MERGE( nest_offl%tind+1, 1, bc_dirichlet_n ), & ! number of elements along x
652	MERGE( nxr-nxl+1, 0, bc_dirichlet_n ), & ! number of elements alogn z
653	MERGE( nest_offl%nzu, 0, bc_dirichlet_n ), & ! number of time steps (2 or 0)
654	MERGE( 2, 0, bc_dirichlet_n ), & ! parallel IO when compiled accordingly
655	.TRUE. )
656	ENDIF
657
658	IF ( air_chemistry ) THEN
659	DO n = 1, UBOUND(nest_offl%var_names_chem_n, 1)
660	IF ( check_existence( nest_offl%var_names, &
661	nest_offl%var_names_chem_n(n) ) ) THEN
662	CALL get_variable( pids_id, &
663	TRIM( nest_offl%var_names_chem_n(n) ), &
664	nest_offl%chem_north(:,:,:,n), &
665	MERGE( nxl+1, 1, bc_dirichlet_n ), &
666	MERGE( nzb+1, 1, bc_dirichlet_n ), &
667	MERGE( nest_offl%tind+1, 1, bc_dirichlet_n ), &
668	MERGE( nxr-nxl+1, 0, bc_dirichlet_n ), &
669	MERGE( nest_offl%nzu, 0, bc_dirichlet_n ), &
670	MERGE( 2, 0, bc_dirichlet_n ), &
671	.TRUE. )
672	nest_offl%chem_from_file_n(n) = .TRUE.
673	ENDIF
674	ENDDO
675	ENDIF
676	!
677	!-- Read data for southern boundary
678	CALL get_variable( pids_id, 'ls_forcing_south_u', & ! array to be read
679	nest_offl%u_south, & ! start index x direction
680	MERGE( nxlu, 1, bc_dirichlet_s ), & ! start index z direction
681	MERGE( nzb+1, 1, bc_dirichlet_s ), & ! start index time dimension
682	MERGE( nest_offl%tind+1, 1, bc_dirichlet_s ), & ! number of elements along x
683	MERGE( nxr-nxlu+1, 0, bc_dirichlet_s ), & ! number of elements alogn z
684	MERGE( nest_offl%nzu, 0, bc_dirichlet_s ), & ! number of time steps (2 or 0)
685	MERGE( 2, 0, bc_dirichlet_s ), & ! parallel IO when compiled accordingly
686	.TRUE. )
687
688	CALL get_variable( pids_id, 'ls_forcing_south_v', & ! array to be read
689	nest_offl%v_south, & ! start index x direction
690	MERGE( nxl+1, 1, bc_dirichlet_s ), & ! start index z direction
691	MERGE( nzb+1, 1, bc_dirichlet_s ), & ! start index time dimension
692	MERGE( nest_offl%tind+1, 1, bc_dirichlet_s ), & ! number of elements along x
693	MERGE( nxr-nxl+1, 0, bc_dirichlet_s ), & ! number of elements alogn z
694	MERGE( nest_offl%nzu, 0, bc_dirichlet_s ), & ! number of time steps (2 or 0)
695	MERGE( 2, 0, bc_dirichlet_s ), & ! parallel IO when compiled accordingly
696	.TRUE. )
697
698	CALL get_variable( pids_id, 'ls_forcing_south_w', & ! array to be read
699	nest_offl%w_south, & ! start index x direction
700	MERGE( nxl+1, 1, bc_dirichlet_s ), & ! start index z direction
701	MERGE( nzb+1, 1, bc_dirichlet_s ), & ! start index time dimension
702	MERGE( nest_offl%tind+1, 1, bc_dirichlet_s ), & ! number of elements along x
703	MERGE( nxr-nxl+1, 0, bc_dirichlet_s ), & ! number of elements alogn z
704	MERGE( nest_offl%nzw, 0, bc_dirichlet_s ), & ! number of time steps (2 or 0)
705	MERGE( 2, 0, bc_dirichlet_s ), & ! parallel IO when compiled accordingly
706	.TRUE. )
707
708	IF ( .NOT. neutral ) THEN
709	CALL get_variable( pids_id, 'ls_forcing_south_pt', & ! array to be read
710	nest_offl%pt_south, & ! start index x direction
711	MERGE( nxl+1, 1, bc_dirichlet_s ), & ! start index z direction
712	MERGE( nzb+1, 1, bc_dirichlet_s ), & ! start index time dimension
713	MERGE( nest_offl%tind+1, 1, bc_dirichlet_s ), & ! number of elements along x
714	MERGE( nxr-nxl+1, 0, bc_dirichlet_s ), & ! number of elements alogn z
715	MERGE( nest_offl%nzu, 0, bc_dirichlet_s ), & ! number of time steps (2 or 0)
716	MERGE( 2, 0, bc_dirichlet_s ), & ! parallel IO when compiled accordingly
717	.TRUE. )
718	ENDIF
719	IF ( humidity ) THEN
720	CALL get_variable( pids_id, 'ls_forcing_south_qv', & ! array to be read
721	nest_offl%q_south, & ! start index x direction
722	MERGE( nxl+1, 1, bc_dirichlet_s ), & ! start index z direction
723	MERGE( nzb+1, 1, bc_dirichlet_s ), & ! start index time dimension
724	MERGE( nest_offl%tind+1, 1, bc_dirichlet_s ), & ! number of elements along x
725	MERGE( nxr-nxl+1, 0, bc_dirichlet_s ), & ! number of elements alogn z
726	MERGE( nest_offl%nzu, 0, bc_dirichlet_s ), & ! number of time steps (2 or 0)
727	MERGE( 2, 0, bc_dirichlet_s ), & ! parallel IO when compiled accordingly
728	.TRUE. )
729	ENDIF
730
731	IF ( air_chemistry ) THEN
732	DO n = 1, UBOUND(nest_offl%var_names_chem_s, 1)
733	IF ( check_existence( nest_offl%var_names, &
734	nest_offl%var_names_chem_s(n) ) ) THEN
735	CALL get_variable( pids_id, &
736	TRIM( nest_offl%var_names_chem_s(n) ), &
737	nest_offl%chem_south(:,:,:,n), &
738	MERGE( nxl+1, 1, bc_dirichlet_s ), &
739	MERGE( nzb+1, 1, bc_dirichlet_s ), &
740	MERGE( nest_offl%tind+1, 1, bc_dirichlet_s ), &
741	MERGE( nxr-nxl+1, 0, bc_dirichlet_s ), &
742	MERGE( nest_offl%nzu, 0, bc_dirichlet_s ), &
743	MERGE( 2, 0, bc_dirichlet_s ), &
744	.TRUE. )
745	nest_offl%chem_from_file_s(n) = .TRUE.
746	ENDIF
747	ENDDO
748	ENDIF
749	!
750	!-- Top boundary
751	CALL get_variable( pids_id, 'ls_forcing_top_u', &
752	nest_offl%u_top(0:1,nys:nyn,nxlu:nxr), &
753	nxlu, nys+1, nest_offl%tind+1, &
754	nxr-nxlu+1, nyn-nys+1, 2, .TRUE. )
755
756	CALL get_variable( pids_id, 'ls_forcing_top_v', &
757	nest_offl%v_top(0:1,nysv:nyn,nxl:nxr), &
758	nxl+1, nysv, nest_offl%tind+1, &
759	nxr-nxl+1, nyn-nysv+1, 2, .TRUE. )
760
761	CALL get_variable( pids_id, 'ls_forcing_top_w', &
762	nest_offl%w_top(0:1,nys:nyn,nxl:nxr), &
763	nxl+1, nys+1, nest_offl%tind+1, &
764	nxr-nxl+1, nyn-nys+1, 2, .TRUE. )
765
766	IF ( .NOT. neutral ) THEN
767	CALL get_variable( pids_id, 'ls_forcing_top_pt', &
768	nest_offl%pt_top(0:1,nys:nyn,nxl:nxr), &
769	nxl+1, nys+1, nest_offl%tind+1, &
770	nxr-nxl+1, nyn-nys+1, 2, .TRUE. )
771	ENDIF
772	IF ( humidity ) THEN
773	CALL get_variable( pids_id, 'ls_forcing_top_qv', &
774	nest_offl%q_top(0:1,nys:nyn,nxl:nxr), &
775	nxl+1, nys+1, nest_offl%tind+1, &
776	nxr-nxl+1, nyn-nys+1, 2, .TRUE. )
777	ENDIF
778
779	IF ( air_chemistry ) THEN
780	DO n = 1, UBOUND(nest_offl%var_names_chem_t, 1)
781	IF ( check_existence( nest_offl%var_names, &
782	nest_offl%var_names_chem_t(n) ) ) THEN
783	CALL get_variable( pids_id, &
784	TRIM( nest_offl%var_names_chem_t(n) ), &
785	nest_offl%chem_top(0:1,nys:nyn,nxl:nxr,n), &
786	nxl+1, nys+1, nest_offl%tind+1, &
787	nxr-nxl+1, nyn-nys+1, 2, .TRUE. )
788	nest_offl%chem_from_file_t(n) = .TRUE.
789	ENDIF
790	ENDDO
791	ENDIF
792
793	!
794	!-- Close input file
795	CALL close_input_file( pids_id )
796	#endif
797	!
798	!-- Set control flag to indicate that boundary data has been initially
799	!-- input.
800	nest_offl%init = .TRUE.
801	!
802	!-- End of CPU measurement
803	CALL cpu_log( log_point_s(86), 'NetCDF input forcing', 'stop' )
804
805	END SUBROUTINE nesting_offl_input
806
807
808	!------------------------------------------------------------------------------!
809	! Description:
810	! ------------
811	!> In this subroutine a constant mass within the model domain is guaranteed.
812	!> Larger-scale models may be based on a compressible equation system, which is
813	!> not consistent with PALMs incompressible equation system. In order to avoid
814	!> a decrease or increase of mass during the simulation, non-divergent flow
815	!> through the lateral and top boundaries is compensated by the vertical wind
816	!> component at the top boundary.
817	!------------------------------------------------------------------------------!
818	SUBROUTINE nesting_offl_mass_conservation
819
820	INTEGER(iwp) :: i !< grid index in x-direction
821	INTEGER(iwp) :: j !< grid index in y-direction
822	INTEGER(iwp) :: k !< grid index in z-direction
823
824	REAL(wp) :: d_area_t !< inverse of the total area of the horizontal model domain
825	REAL(wp) :: w_correct !< vertical velocity increment required to compensate non-divergent flow through the boundaries
826	REAL(wp), DIMENSION(1:3) :: volume_flow_l !< local volume flow
827
828
829	IF ( debug_output_timestep ) CALL debug_message( 'nesting_offl_mass_conservation', 'start' )
830
831	CALL cpu_log( log_point(58), 'offline nesting', 'start' )
832
833	volume_flow = 0.0_wp
834	volume_flow_l = 0.0_wp
835
836	d_area_t = 1.0_wp / ( ( nx + 1 ) * dx * ( ny + 1 ) * dy )
837
838	IF ( bc_dirichlet_l ) THEN
839	i = nxl
840	DO j = nys, nyn
841	DO k = nzb+1, nzt
842	volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k) * dy &
843	* MERGE( 1.0_wp, 0.0_wp, &
844	BTEST( wall_flags_0(k,j,i), 1 ) )
845	ENDDO
846	ENDDO
847	ENDIF
848	IF ( bc_dirichlet_r ) THEN
849	i = nxr+1
850	DO j = nys, nyn
851	DO k = nzb+1, nzt
852	volume_flow_l(1) = volume_flow_l(1) - u(k,j,i) * dzw(k) * dy &
853	* MERGE( 1.0_wp, 0.0_wp, &
854	BTEST( wall_flags_0(k,j,i), 1 ) )
855	ENDDO
856	ENDDO
857	ENDIF
858	IF ( bc_dirichlet_s ) THEN
859	j = nys
860	DO i = nxl, nxr
861	DO k = nzb+1, nzt
862	volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k) * dx &
863	* MERGE( 1.0_wp, 0.0_wp, &
864	BTEST( wall_flags_0(k,j,i), 2 ) )
865	ENDDO
866	ENDDO
867	ENDIF
868	IF ( bc_dirichlet_n ) THEN
869	j = nyn+1
870	DO i = nxl, nxr
871	DO k = nzb+1, nzt
872	volume_flow_l(2) = volume_flow_l(2) - v(k,j,i) * dzw(k) * dx &
873	* MERGE( 1.0_wp, 0.0_wp, &
874	BTEST( wall_flags_0(k,j,i), 2 ) )
875	ENDDO
876	ENDDO
877	ENDIF
878	!
879	!-- Top boundary
880	k = nzt
881	DO i = nxl, nxr
882	DO j = nys, nyn
883	volume_flow_l(3) = volume_flow_l(3) - w(k,j,i) * dx * dy
884	ENDDO
885	ENDDO
886
887	#if defined( __parallel )
888	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
889	CALL MPI_ALLREDUCE( volume_flow_l, volume_flow, 3, MPI_REAL, MPI_SUM, &
890	comm2d, ierr )
891	#else
892	volume_flow = volume_flow_l
893	#endif
894
895	w_correct = SUM( volume_flow ) * d_area_t
896
897	DO i = nxl, nxr
898	DO j = nys, nyn
899	DO k = nzt, nzt + 1
900	w(k,j,i) = w(k,j,i) + w_correct
901	ENDDO
902	ENDDO
903	ENDDO
904
905	CALL cpu_log( log_point(58), 'offline nesting', 'stop' )
906
907	IF ( debug_output_timestep ) CALL debug_message( 'nesting_offl_mass_conservation', 'end' )
908
909	END SUBROUTINE nesting_offl_mass_conservation
910
911
912	!------------------------------------------------------------------------------!
913	! Description:
914	! ------------
915	!> Set the lateral and top boundary conditions in case the PALM domain is
916	!> nested offline in a mesoscale model. Further, average boundary data and
917	!> determine mean profiles, further used for correct damping in the sponge
918	!> layer.
919	!------------------------------------------------------------------------------!
920	SUBROUTINE nesting_offl_bc
921
922	INTEGER(iwp) :: i !< running index x-direction
923	INTEGER(iwp) :: j !< running index y-direction
924	INTEGER(iwp) :: k !< running index z-direction
925	INTEGER(iwp) :: n !< running index for chemical species
926
927	REAL(wp), DIMENSION(nzb:nzt+1) :: pt_ref !< reference profile for potential temperature
928	REAL(wp), DIMENSION(nzb:nzt+1) :: pt_ref_l !< reference profile for potential temperature on subdomain
929	REAL(wp), DIMENSION(nzb:nzt+1) :: q_ref !< reference profile for mixing ratio
930	REAL(wp), DIMENSION(nzb:nzt+1) :: q_ref_l !< reference profile for mixing ratio on subdomain
931	REAL(wp), DIMENSION(nzb:nzt+1) :: u_ref !< reference profile for u-component
932	REAL(wp), DIMENSION(nzb:nzt+1) :: u_ref_l !< reference profile for u-component on subdomain
933	REAL(wp), DIMENSION(nzb:nzt+1) :: v_ref !< reference profile for v-component
934	REAL(wp), DIMENSION(nzb:nzt+1) :: v_ref_l !< reference profile for v-component on subdomain
935
936	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ref_chem !< reference profile for chemical species
937	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ref_chem_l !< reference profile for chemical species on subdomain
938
939	IF ( debug_output_timestep ) CALL debug_message( 'nesting_offl_bc', 'start' )
940
941	CALL cpu_log( log_point(58), 'offline nesting', 'start' )
942	!
943	!-- Initialize mean profiles, derived from boundary data, to zero
944	pt_ref = 0.0_wp
945	q_ref = 0.0_wp
946	u_ref = 0.0_wp
947	v_ref = 0.0_wp
948
949	pt_ref_l = 0.0_wp
950	q_ref_l = 0.0_wp
951	u_ref_l = 0.0_wp
952	v_ref_l = 0.0_wp
953	!
954	!-- If required, allocate temporary arrays to compute chemistry mean profiles
955	IF ( air_chemistry ) THEN
956	ALLOCATE( ref_chem(nzb:nzt+1,1:UBOUND( chem_species, 1 ) ) )
957	ALLOCATE( ref_chem_l(nzb:nzt+1,1:UBOUND( chem_species, 1 ) ) )
958	ref_chem = 0.0_wp
959	ref_chem_l = 0.0_wp
960	ENDIF
961	!
962	!-- Set boundary conditions of u-, v-, w-component, as well as q, and pt.
963	!-- Note, boundary values at the left boundary: i=-1 (v,w,pt,q) and
964	!-- i=0 (u), at the right boundary: i=nxr+1 (all), at the south boundary:
965	!-- j=-1 (u,w,pt,q) and j=0 (v), at the north boundary: j=nyn+1 (all).
966	!-- Please note, at the left (for u) and south (for v) boundary, values
967	!-- for u and v are set also at i/j=-1, since these values are used in
968	!-- boundary_conditions() to restore prognostic values.
969	!-- Further, sum up data to calculate mean profiles from boundary data,
970	!-- used for Rayleigh damping.
971	IF ( bc_dirichlet_l ) THEN
972
973	DO j = nys, nyn
974	DO k = nzb+1, nzt
975	u(k,j,0) = interpolate_in_time( nest_offl%u_left(0,k,j), &
976	nest_offl%u_left(1,k,j), &
977	fac_dt ) * &
978	MERGE( 1.0_wp, 0.0_wp, &
979	BTEST( wall_flags_0(k,j,0), 1 ) )
980	u(k,j,-1) = u(k,j,0)
981	ENDDO
982	u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,j,0)
983	ENDDO
984
985	DO j = nys, nyn
986	DO k = nzb+1, nzt-1
987	w(k,j,-1) = interpolate_in_time( nest_offl%w_left(0,k,j), &
988	nest_offl%w_left(1,k,j), &
989	fac_dt ) * &
990	MERGE( 1.0_wp, 0.0_wp, &
991	BTEST( wall_flags_0(k,j,-1), 3 ) )
992	ENDDO
993	w(nzt,j,-1) = w(nzt-1,j,-1)
994	ENDDO
995
996	DO j = nysv, nyn
997	DO k = nzb+1, nzt
998	v(k,j,-1) = interpolate_in_time( nest_offl%v_left(0,k,j), &
999	nest_offl%v_left(1,k,j), &
1000	fac_dt ) * &
1001	MERGE( 1.0_wp, 0.0_wp, &
1002	BTEST( wall_flags_0(k,j,-1), 2 ) )
1003	ENDDO
1004	v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,j,-1)
1005	ENDDO
1006
1007	IF ( .NOT. neutral ) THEN
1008	DO j = nys, nyn
1009	DO k = nzb+1, nzt
1010	pt(k,j,-1) = interpolate_in_time( nest_offl%pt_left(0,k,j), &
1011	nest_offl%pt_left(1,k,j), &
1012	fac_dt )
1013
1014	ENDDO
1015	pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,j,-1)
1016	ENDDO
1017	ENDIF
1018
1019	IF ( humidity ) THEN
1020	DO j = nys, nyn
1021	DO k = nzb+1, nzt
1022	q(k,j,-1) = interpolate_in_time( nest_offl%q_left(0,k,j), &
1023	nest_offl%q_left(1,k,j), &
1024	fac_dt )
1025
1026	ENDDO
1027	q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,j,-1)
1028	ENDDO
1029	ENDIF
1030
1031	IF ( air_chemistry ) THEN
1032	DO n = 1, UBOUND( chem_species, 1 )
1033	IF ( nest_offl%chem_from_file_l(n) ) THEN
1034	DO j = nys, nyn
1035	DO k = nzb+1, nzt
1036	chem_species(n)%conc(k,j,-1) = interpolate_in_time( &
1037	nest_offl%chem_left(0,k,j,n),&
1038	nest_offl%chem_left(1,k,j,n),&
1039	fac_dt )
1040	ENDDO
1041	ref_chem_l(nzb+1:nzt,n) = ref_chem_l(nzb+1:nzt,n) &
1042	+ chem_species(n)%conc(nzb+1:nzt,j,-1)
1043	ENDDO
1044	ENDIF
1045	ENDDO
1046	ENDIF
1047
1048	ENDIF
1049
1050	IF ( bc_dirichlet_r ) THEN
1051
1052	DO j = nys, nyn
1053	DO k = nzb+1, nzt
1054	u(k,j,nxr+1) = interpolate_in_time( nest_offl%u_right(0,k,j), &
1055	nest_offl%u_right(1,k,j), &
1056	fac_dt ) * &
1057	MERGE( 1.0_wp, 0.0_wp, &
1058	BTEST( wall_flags_0(k,j,nxr+1), 1 ) )
1059	ENDDO
1060	u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,j,nxr+1)
1061	ENDDO
1062	DO j = nys, nyn
1063	DO k = nzb+1, nzt-1
1064	w(k,j,nxr+1) = interpolate_in_time( nest_offl%w_right(0,k,j), &
1065	nest_offl%w_right(1,k,j), &
1066	fac_dt ) * &
1067	MERGE( 1.0_wp, 0.0_wp, &
1068	BTEST( wall_flags_0(k,j,nxr+1), 3 ) )
1069	ENDDO
1070	w(nzt,j,nxr+1) = w(nzt-1,j,nxr+1)
1071	ENDDO
1072
1073	DO j = nysv, nyn
1074	DO k = nzb+1, nzt
1075	v(k,j,nxr+1) = interpolate_in_time( nest_offl%v_right(0,k,j), &
1076	nest_offl%v_right(1,k,j), &
1077	fac_dt ) * &
1078	MERGE( 1.0_wp, 0.0_wp, &
1079	BTEST( wall_flags_0(k,j,nxr+1), 2 ) )
1080	ENDDO
1081	v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,j,nxr+1)
1082	ENDDO
1083
1084	IF ( .NOT. neutral ) THEN
1085	DO j = nys, nyn
1086	DO k = nzb+1, nzt
1087	pt(k,j,nxr+1) = interpolate_in_time( &
1088	nest_offl%pt_right(0,k,j), &
1089	nest_offl%pt_right(1,k,j), &
1090	fac_dt )
1091	ENDDO
1092	pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,j,nxr+1)
1093	ENDDO
1094	ENDIF
1095
1096	IF ( humidity ) THEN
1097	DO j = nys, nyn
1098	DO k = nzb+1, nzt
1099	q(k,j,nxr+1) = interpolate_in_time( &
1100	nest_offl%q_right(0,k,j), &
1101	nest_offl%q_right(1,k,j), &
1102	fac_dt )
1103
1104	ENDDO
1105	q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,j,nxr+1)
1106	ENDDO
1107	ENDIF
1108
1109	IF ( air_chemistry ) THEN
1110	DO n = 1, UBOUND( chem_species, 1 )
1111	IF ( nest_offl%chem_from_file_r(n) ) THEN
1112	DO j = nys, nyn
1113	DO k = nzb+1, nzt
1114	chem_species(n)%conc(k,j,nxr+1) = interpolate_in_time(&
1115	nest_offl%chem_right(0,k,j,n),&
1116	nest_offl%chem_right(1,k,j,n),&
1117	fac_dt )
1118	ENDDO
1119	ref_chem_l(nzb+1:nzt,n) = ref_chem_l(nzb+1:nzt,n) &
1120	+ chem_species(n)%conc(nzb+1:nzt,j,nxr+1)
1121	ENDDO
1122	ENDIF
1123	ENDDO
1124	ENDIF
1125
1126	ENDIF
1127
1128	IF ( bc_dirichlet_s ) THEN
1129
1130	DO i = nxl, nxr
1131	DO k = nzb+1, nzt
1132	v(k,0,i) = interpolate_in_time( nest_offl%v_south(0,k,i), &
1133	nest_offl%v_south(1,k,i), &
1134	fac_dt ) * &
1135	MERGE( 1.0_wp, 0.0_wp, &
1136	BTEST( wall_flags_0(k,0,i), 2 ) )
1137	v(k,-1,i) = v(k,0,i)
1138	ENDDO
1139	v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,0,i)
1140	ENDDO
1141
1142	DO i = nxl, nxr
1143	DO k = nzb+1, nzt-1
1144	w(k,-1,i) = interpolate_in_time( nest_offl%w_south(0,k,i), &
1145	nest_offl%w_south(1,k,i), &
1146	fac_dt ) * &
1147	MERGE( 1.0_wp, 0.0_wp, &
1148	BTEST( wall_flags_0(k,-1,i), 3 ) )
1149	ENDDO
1150	w(nzt,-1,i) = w(nzt-1,-1,i)
1151	ENDDO
1152
1153	DO i = nxlu, nxr
1154	DO k = nzb+1, nzt
1155	u(k,-1,i) = interpolate_in_time( nest_offl%u_south(0,k,i), &
1156	nest_offl%u_south(1,k,i), &
1157	fac_dt ) * &
1158	MERGE( 1.0_wp, 0.0_wp, &
1159	BTEST( wall_flags_0(k,-1,i), 1 ) )
1160	ENDDO
1161	u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,-1,i)
1162	ENDDO
1163
1164	IF ( .NOT. neutral ) THEN
1165	DO i = nxl, nxr
1166	DO k = nzb+1, nzt
1167	pt(k,-1,i) = interpolate_in_time( &
1168	nest_offl%pt_south(0,k,i), &
1169	nest_offl%pt_south(1,k,i), &
1170	fac_dt )
1171
1172	ENDDO
1173	pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,-1,i)
1174	ENDDO
1175	ENDIF
1176
1177	IF ( humidity ) THEN
1178	DO i = nxl, nxr
1179	DO k = nzb+1, nzt
1180	q(k,-1,i) = interpolate_in_time( &
1181	nest_offl%q_south(0,k,i), &
1182	nest_offl%q_south(1,k,i), &
1183	fac_dt )
1184
1185	ENDDO
1186	q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,-1,i)
1187	ENDDO
1188	ENDIF
1189
1190	IF ( air_chemistry ) THEN
1191	DO n = 1, UBOUND( chem_species, 1 )
1192	IF ( nest_offl%chem_from_file_s(n) ) THEN
1193	DO i = nxl, nxr
1194	DO k = nzb+1, nzt
1195	chem_species(n)%conc(k,-1,i) = interpolate_in_time( &
1196	nest_offl%chem_south(0,k,i,n),&
1197	nest_offl%chem_south(1,k,i,n),&
1198	fac_dt )
1199	ENDDO
1200	ref_chem_l(nzb+1:nzt,n) = ref_chem_l(nzb+1:nzt,n) &
1201	+ chem_species(n)%conc(nzb+1:nzt,-1,i)
1202	ENDDO
1203	ENDIF
1204	ENDDO
1205	ENDIF
1206
1207	ENDIF
1208
1209	IF ( bc_dirichlet_n ) THEN
1210
1211	DO i = nxl, nxr
1212	DO k = nzb+1, nzt
1213	v(k,nyn+1,i) = interpolate_in_time( nest_offl%v_north(0,k,i), &
1214	nest_offl%v_north(1,k,i), &
1215	fac_dt ) * &
1216	MERGE( 1.0_wp, 0.0_wp, &
1217	BTEST( wall_flags_0(k,nyn+1,i), 2 ) )
1218	ENDDO
1219	v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,nyn+1,i)
1220	ENDDO
1221	DO i = nxl, nxr
1222	DO k = nzb+1, nzt-1
1223	w(k,nyn+1,i) = interpolate_in_time( nest_offl%w_north(0,k,i), &
1224	nest_offl%w_north(1,k,i), &
1225	fac_dt ) * &
1226	MERGE( 1.0_wp, 0.0_wp, &
1227	BTEST( wall_flags_0(k,nyn+1,i), 3 ) )
1228	ENDDO
1229	w(nzt,nyn+1,i) = w(nzt-1,nyn+1,i)
1230	ENDDO
1231
1232	DO i = nxlu, nxr
1233	DO k = nzb+1, nzt
1234	u(k,nyn+1,i) = interpolate_in_time( nest_offl%u_north(0,k,i), &
1235	nest_offl%u_north(1,k,i), &
1236	fac_dt ) * &
1237	MERGE( 1.0_wp, 0.0_wp, &
1238	BTEST( wall_flags_0(k,nyn+1,i), 1 ) )
1239
1240	ENDDO
1241	u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,nyn+1,i)
1242	ENDDO
1243
1244	IF ( .NOT. neutral ) THEN
1245	DO i = nxl, nxr
1246	DO k = nzb+1, nzt
1247	pt(k,nyn+1,i) = interpolate_in_time( &
1248	nest_offl%pt_north(0,k,i), &
1249	nest_offl%pt_north(1,k,i), &
1250	fac_dt )
1251
1252	ENDDO
1253	pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,nyn+1,i)
1254	ENDDO
1255	ENDIF
1256
1257	IF ( humidity ) THEN
1258	DO i = nxl, nxr
1259	DO k = nzb+1, nzt
1260	q(k,nyn+1,i) = interpolate_in_time( &
1261	nest_offl%q_north(0,k,i), &
1262	nest_offl%q_north(1,k,i), &
1263	fac_dt )
1264
1265	ENDDO
1266	q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,nyn+1,i)
1267	ENDDO
1268	ENDIF
1269
1270	IF ( air_chemistry ) THEN
1271	DO n = 1, UBOUND( chem_species, 1 )
1272	IF ( nest_offl%chem_from_file_n(n) ) THEN
1273	DO i = nxl, nxr
1274	DO k = nzb+1, nzt
1275	chem_species(n)%conc(k,nyn+1,i) = interpolate_in_time(&
1276	nest_offl%chem_north(0,k,i,n),&
1277	nest_offl%chem_north(1,k,i,n),&
1278	fac_dt )
1279	ENDDO
1280	ref_chem_l(nzb+1:nzt,n) = ref_chem_l(nzb+1:nzt,n) &
1281	+ chem_species(n)%conc(nzb+1:nzt,nyn+1,i)
1282	ENDDO
1283	ENDIF
1284	ENDDO
1285	ENDIF
1286
1287	ENDIF
1288	!
1289	!-- Top boundary
1290	DO i = nxlu, nxr
1291	DO j = nys, nyn
1292	u(nzt+1,j,i) = interpolate_in_time( nest_offl%u_top(0,j,i), &
1293	nest_offl%u_top(1,j,i), &
1294	fac_dt ) * &
1295	MERGE( 1.0_wp, 0.0_wp, &
1296	BTEST( wall_flags_0(nzt+1,j,i), 1 ) )
1297	u_ref_l(nzt+1) = u_ref_l(nzt+1) + u(nzt+1,j,i)
1298	ENDDO
1299	ENDDO
1300	!
1301	!-- For left boundary set boundary condition for u-component also at top
1302	!-- grid point.
1303	!-- Note, this has no effect on the numeric solution, only for data output.
1304	IF ( bc_dirichlet_l ) u(nzt+1,:,nxl) = u(nzt+1,:,nxlu)
1305
1306	DO i = nxl, nxr
1307	DO j = nysv, nyn
1308	v(nzt+1,j,i) = interpolate_in_time( nest_offl%v_top(0,j,i), &
1309	nest_offl%v_top(1,j,i), &
1310	fac_dt ) * &
1311	MERGE( 1.0_wp, 0.0_wp, &
1312	BTEST( wall_flags_0(nzt+1,j,i), 2 ) )
1313	v_ref_l(nzt+1) = v_ref_l(nzt+1) + v(nzt+1,j,i)
1314	ENDDO
1315	ENDDO
1316	!
1317	!-- For south boundary set boundary condition for v-component also at top
1318	!-- grid point.
1319	!-- Note, this has no effect on the numeric solution, only for data output.
1320	IF ( bc_dirichlet_s ) v(nzt+1,nys,:) = v(nzt+1,nysv,:)
1321
1322	DO i = nxl, nxr
1323	DO j = nys, nyn
1324	w(nzt,j,i) = interpolate_in_time( nest_offl%w_top(0,j,i), &
1325	nest_offl%w_top(1,j,i), &
1326	fac_dt ) * &
1327	MERGE( 1.0_wp, 0.0_wp, &
1328	BTEST( wall_flags_0(nzt,j,i), 3 ) )
1329	w(nzt+1,j,i) = w(nzt,j,i)
1330	ENDDO
1331	ENDDO
1332
1333
1334	IF ( .NOT. neutral ) THEN
1335	DO i = nxl, nxr
1336	DO j = nys, nyn
1337	pt(nzt+1,j,i) = interpolate_in_time( nest_offl%pt_top(0,j,i), &
1338	nest_offl%pt_top(1,j,i), &
1339	fac_dt )
1340	pt_ref_l(nzt+1) = pt_ref_l(nzt+1) + pt(nzt+1,j,i)
1341	ENDDO
1342	ENDDO
1343	ENDIF
1344
1345	IF ( humidity ) THEN
1346	DO i = nxl, nxr
1347	DO j = nys, nyn
1348	q(nzt+1,j,i) = interpolate_in_time( nest_offl%q_top(0,j,i), &
1349	nest_offl%q_top(1,j,i), &
1350	fac_dt )
1351	q_ref_l(nzt+1) = q_ref_l(nzt+1) + q(nzt+1,j,i)
1352	ENDDO
1353	ENDDO
1354	ENDIF
1355
1356	IF ( air_chemistry ) THEN
1357	DO n = 1, UBOUND( chem_species, 1 )
1358	IF ( nest_offl%chem_from_file_t(n) ) THEN
1359	DO i = nxl, nxr
1360	DO j = nys, nyn
1361	chem_species(n)%conc(nzt+1,j,i) = interpolate_in_time( &
1362	nest_offl%chem_top(0,j,i,n), &
1363	nest_offl%chem_top(1,j,i,n), &
1364	fac_dt )
1365	ref_chem_l(nzt+1,n) = ref_chem_l(nzt+1,n) + &
1366	chem_species(n)%conc(nzt+1,j,i)
1367	ENDDO
1368	ENDDO
1369	ENDIF
1370	ENDDO
1371	ENDIF
1372	!
1373	!-- Moreover, set Neumann boundary condition for subgrid-scale TKE,
1374	!-- passive scalar, dissipation, and chemical species if required
1375	IF ( rans_mode .AND. rans_tke_e ) THEN
1376	IF ( bc_dirichlet_l ) diss(:,:,nxl-1) = diss(:,:,nxl)
1377	IF ( bc_dirichlet_r ) diss(:,:,nxr+1) = diss(:,:,nxr)
1378	IF ( bc_dirichlet_s ) diss(:,nys-1,:) = diss(:,nys,:)
1379	IF ( bc_dirichlet_n ) diss(:,nyn+1,:) = diss(:,nyn,:)
1380	ENDIF
1381	! IF ( .NOT. constant_diffusion ) THEN
1382	! IF ( bc_dirichlet_l ) e(:,:,nxl-1) = e(:,:,nxl)
1383	! IF ( bc_dirichlet_r ) e(:,:,nxr+1) = e(:,:,nxr)
1384	! IF ( bc_dirichlet_s ) e(:,nys-1,:) = e(:,nys,:)
1385	! IF ( bc_dirichlet_n ) e(:,nyn+1,:) = e(:,nyn,:)
1386	! e(nzt+1,:,:) = e(nzt,:,:)
1387	! ENDIF
1388	! IF ( passive_scalar ) THEN
1389	! IF ( bc_dirichlet_l ) s(:,:,nxl-1) = s(:,:,nxl)
1390	! IF ( bc_dirichlet_r ) s(:,:,nxr+1) = s(:,:,nxr)
1391	! IF ( bc_dirichlet_s ) s(:,nys-1,:) = s(:,nys,:)
1392	! IF ( bc_dirichlet_n ) s(:,nyn+1,:) = s(:,nyn,:)
1393	! ENDIF
1394
1395	CALL exchange_horiz( u, nbgp )
1396	CALL exchange_horiz( v, nbgp )
1397	CALL exchange_horiz( w, nbgp )
1398	IF ( .NOT. neutral ) CALL exchange_horiz( pt, nbgp )
1399	IF ( humidity ) CALL exchange_horiz( q, nbgp )
1400	IF ( air_chemistry ) THEN
1401	DO n = 1, UBOUND( chem_species, 1 )
1402	!
1403	!-- Do local exchange only when necessary, i.e. when data is coming
1404	!-- from dynamic file.
1405	IF ( nest_offl%chem_from_file_t(n) ) &
1406	CALL exchange_horiz( chem_species(n)%conc, nbgp )
1407	ENDDO
1408	ENDIF
1409	!
1410	!-- Set top boundary condition at all horizontal grid points, also at the
1411	!-- lateral boundary grid points.
1412	w(nzt+1,:,:) = w(nzt,:,:)
1413	!
1414	!-- In case of Rayleigh damping, where the profiles u_init, v_init
1415	!-- q_init and pt_init are still used, update these profiles from the
1416	!-- averaged boundary data.
1417	!-- But first, average these data.
1418	#if defined( __parallel )
1419	CALL MPI_ALLREDUCE( u_ref_l, u_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM, &
1420	comm2d, ierr )
1421	CALL MPI_ALLREDUCE( v_ref_l, v_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM, &
1422	comm2d, ierr )
1423	IF ( humidity ) THEN
1424	CALL MPI_ALLREDUCE( q_ref_l, q_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM, &
1425	comm2d, ierr )
1426	ENDIF
1427	IF ( .NOT. neutral ) THEN
1428	CALL MPI_ALLREDUCE( pt_ref_l, pt_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM,&
1429	comm2d, ierr )
1430	ENDIF
1431	IF ( air_chemistry ) THEN
1432	CALL MPI_ALLREDUCE( ref_chem_l, ref_chem, &
1433	( nzt+1-nzb+1 ) * SIZE( ref_chem(nzb,:) ), &
1434	MPI_REAL, MPI_SUM, comm2d, ierr )
1435	ENDIF
1436	#else
1437	u_ref = u_ref_l
1438	v_ref = v_ref_l
1439	IF ( humidity ) q_ref = q_ref_l
1440	IF ( .NOT. neutral ) pt_ref = pt_ref_l
1441	IF ( air_chemistry ) ref_chem = ref_chem_l
1442	#endif
1443	!
1444	!-- Average data. Note, reference profiles up to nzt are derived from lateral
1445	!-- boundaries, at the model top it is derived from the top boundary. Thus,
1446	!-- number of input data is different from nzb:nzt compared to nzt+1.
1447	!-- Derived from lateral boundaries.
1448	u_ref(nzb:nzt) = u_ref(nzb:nzt) / REAL( 2.0_wp * ( ny + 1 + nx ), &
1449	KIND = wp )
1450	v_ref(nzb:nzt) = v_ref(nzb:nzt) / REAL( 2.0_wp * ( ny + nx + 1 ), &
1451	KIND = wp )
1452	IF ( humidity ) &
1453	q_ref(nzb:nzt) = q_ref(nzb:nzt) / REAL( 2.0_wp * &
1454	( ny + 1 + nx + 1 ), &
1455	KIND = wp )
1456	IF ( .NOT. neutral ) &
1457	pt_ref(nzb:nzt) = pt_ref(nzb:nzt) / REAL( 2.0_wp * &
1458	( ny + 1 + nx + 1 ), &
1459	KIND = wp )
1460	IF ( air_chemistry ) &
1461	ref_chem(nzb:nzt,:) = ref_chem(nzb:nzt,:) / REAL( 2.0_wp * &
1462	( ny + 1 + nx + 1 ), &
1463	KIND = wp )
1464	!
1465	!-- Derived from top boundary.
1466	u_ref(nzt+1) = u_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx ), KIND = wp )
1467	v_ref(nzt+1) = v_ref(nzt+1) / REAL( ( ny ) * ( nx + 1 ), KIND = wp )
1468	IF ( humidity ) &
1469	q_ref(nzt+1) = q_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx + 1 ), &
1470	KIND = wp )
1471	IF ( .NOT. neutral ) &
1472	pt_ref(nzt+1) = pt_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx + 1 ), &
1473	KIND = wp )
1474	IF ( air_chemistry ) &
1475	ref_chem(nzt+1,:) = ref_chem(nzt+1,:) / &
1476	REAL( ( ny + 1 ) * ( nx + 1 ),KIND = wp )
1477	!
1478	!-- Write onto init profiles, which are used for damping. Also set lower
1479	!-- boundary condition for scalars (not required for u and v as these are
1480	!-- zero at k=nzb.
1481	u_init = u_ref
1482	v_init = v_ref
1483	IF ( humidity ) THEN
1484	q_init = q_ref
1485	q_init(nzb) = q_init(nzb+1)
1486	ENDIF
1487	IF ( .NOT. neutral ) THEN
1488	pt_init = pt_ref
1489	pt_init(nzb) = pt_init(nzb+1)
1490	ENDIF
1491
1492	IF ( air_chemistry ) THEN
1493	DO n = 1, UBOUND( chem_species, 1 )
1494	IF ( nest_offl%chem_from_file_t(n) ) THEN
1495	chem_species(n)%conc_pr_init(:) = ref_chem(:,n)
1496	chem_species(n)%conc_pr_init(nzb) = &
1497	chem_species(n)%conc_pr_init(nzb+1)
1498	ENDIF
1499	ENDDO
1500	ENDIF
1501
1502	DEALLOCATE( ref_chem )
1503	DEALLOCATE( ref_chem_l )
1504	!
1505	!-- Further, adjust Rayleigh damping height in case of time-changing conditions.
1506	!-- Therefore, calculate boundary-layer depth first.
1507	CALL nesting_offl_calc_zi
1508	CALL adjust_sponge_layer
1509
1510	CALL cpu_log( log_point(58), 'offline nesting', 'stop' )
1511
1512	IF ( debug_output_timestep ) CALL debug_message( 'nesting_offl_bc', 'end' )
1513
1514
1515	END SUBROUTINE nesting_offl_bc
1516
1517	!------------------------------------------------------------------------------!
1518	! Description:
1519	!------------------------------------------------------------------------------!
1520	!> Update of the geostrophic wind components.
1521	!> @todo: update geostrophic wind also in the child domains (should be done
1522	!> in the nesting.
1523	!------------------------------------------------------------------------------!
1524	SUBROUTINE nesting_offl_geostrophic_wind
1525
1526	INTEGER(iwp) :: k
1527	!
1528	!-- Update geostrophic wind components from dynamic input file.
1529	DO k = nzb+1, nzt
1530	ug(k) = interpolate_in_time( nest_offl%ug(0,k), nest_offl%ug(1,k), &
1531	fac_dt )
1532	vg(k) = interpolate_in_time( nest_offl%vg(0,k), nest_offl%vg(1,k), &
1533	fac_dt )
1534	ENDDO
1535	ug(nzt+1) = ug(nzt)
1536	vg(nzt+1) = vg(nzt)
1537
1538	END SUBROUTINE nesting_offl_geostrophic_wind
1539
1540	!------------------------------------------------------------------------------!
1541	! Description:
1542	!------------------------------------------------------------------------------!
1543	!> Determine the interpolation constant for time interpolation. The
1544	!> calculation is separated from the nesting_offl_bc and
1545	!> nesting_offl_geostrophic_wind in order to be independent on the order
1546	!> of calls.
1547	!------------------------------------------------------------------------------!
1548	SUBROUTINE nesting_offl_interpolation_factor
1549	!
1550	!-- Determine interpolation factor and limit it to 1. This is because
1551	!-- t+dt can slightly exceed time(tind_p) before boundary data is updated
1552	!-- again.
1553	fac_dt = ( time_utc_init + time_since_reference_point &
1554	- nest_offl%time(nest_offl%tind) + dt_3d ) / &
1555	( nest_offl%time(nest_offl%tind_p) - nest_offl%time(nest_offl%tind) )
1556
1557	fac_dt = MIN( 1.0_wp, fac_dt )
1558
1559	END SUBROUTINE nesting_offl_interpolation_factor
1560
1561	!------------------------------------------------------------------------------!
1562	! Description:
1563	!------------------------------------------------------------------------------!
1564	!> Calculates the boundary-layer depth from the boundary data, according to
1565	!> bulk-Richardson criterion.
1566	!------------------------------------------------------------------------------!
1567	SUBROUTINE nesting_offl_calc_zi
1568
1569	INTEGER(iwp) :: i !< loop index in x-direction
1570	INTEGER(iwp) :: j !< loop index in y-direction
1571	INTEGER(iwp) :: k !< loop index in z-direction
1572	INTEGER(iwp) :: k_max_loc !< index of maximum wind speed along z-direction
1573	INTEGER(iwp) :: k_surface !< topography top index in z-direction
1574	INTEGER(iwp) :: num_boundary_gp_non_cyclic !< number of non-cyclic boundaries, used for averaging ABL depth
1575	INTEGER(iwp) :: num_boundary_gp_non_cyclic_l !< number of non-cyclic boundaries, used for averaging ABL depth
1576
1577	REAL(wp) :: ri_bulk !< bulk Richardson number
1578	REAL(wp) :: ri_bulk_crit = 0.25_wp !< critical bulk Richardson number
1579	REAL(wp) :: topo_max !< maximum topography level in model domain
1580	REAL(wp) :: topo_max_l !< maximum topography level in subdomain
1581	REAL(wp) :: vpt_surface !< near-surface virtual potential temperature
1582	REAL(wp) :: zi_l !< mean boundary-layer depth on subdomain
1583	REAL(wp) :: zi_local !< local boundary-layer depth
1584
1585	REAL(wp), DIMENSION(nzb:nzt+1) :: vpt_col !< vertical profile of virtual potential temperature at (j,i)-grid point
1586	REAL(wp), DIMENSION(nzb:nzt+1) :: uv_abs !< vertical profile of horizontal wind speed at (j,i)-grid point
1587
1588
1589	!
1590	!-- Calculate mean boundary-layer height from boundary data.
1591	!-- Start with the left and right boundaries.
1592	zi_l = 0.0_wp
1593	num_boundary_gp_non_cyclic_l = 0
1594	IF ( bc_dirichlet_l .OR. bc_dirichlet_r ) THEN
1595	!
1596	!-- Sum-up and store number of boundary grid points used for averaging
1597	!-- ABL depth
1598	num_boundary_gp_non_cyclic_l = num_boundary_gp_non_cyclic_l + &
1599	nxr - nxl + 1
1600	!
1601	!-- Determine index along x. Please note, index indicates boundary
1602	!-- grid point for scalars.
1603	i = MERGE( -1, nxr + 1, bc_dirichlet_l )
1604
1605	DO j = nys, nyn
1606	!
1607	!-- Determine topography top index at current (j,i) index
1608	k_surface = topo_top_ind(j,i,0)
1609	!
1610	!-- Pre-compute surface virtual temperature. Therefore, use 2nd
1611	!-- prognostic level according to Heinze et al. (2017).
1612	IF ( humidity ) THEN
1613	vpt_surface = pt(k_surface+2,j,i) * &
1614	( 1.0_wp + 0.61_wp * q(k_surface+2,j,i) )
1615	vpt_col = pt(:,j,i) * ( 1.0_wp + 0.61_wp * q(:,j,i) )
1616	ELSE
1617	vpt_surface = pt(k_surface+2,j,i)
1618	vpt_col = pt(:,j,i)
1619	ENDIF
1620	!
1621	!-- Calculate local boundary layer height from bulk Richardson number,
1622	!-- i.e. the height where the bulk Richardson number exceeds its
1623	!-- critical value of 0.25 (according to Heinze et al., 2017).
1624	!-- Note, no interpolation of u- and v-component is made, as both
1625	!-- are mainly mean inflow profiles with very small spatial variation.
1626	!-- Add a safety factor in case the velocity term becomes zero. This
1627	!-- may happen if overhanging 3D structures are directly located at
1628	!-- the boundary, where velocity inside the building is zero
1629	!-- (k_surface is the index of the lowest upward-facing surface).
1630	uv_abs(:) = SQRT( MERGE( u(:,j,i+1), u(:,j,i), &
1631	bc_dirichlet_l )**2 + &
1632	v(:,j,i)**2 )
1633	!
1634	!-- Determine index of the maximum wind speed
1635	k_max_loc = MAXLOC( uv_abs(:), DIM = 1 ) - 1
1636
1637	zi_local = 0.0_wp
1638	DO k = k_surface+1, nzt
1639	ri_bulk = zu(k) * g / vpt_surface * &
1640	( vpt_col(k) - vpt_surface ) / &
1641	( uv_abs(k) + 1E-5_wp )
1642	!
1643	!-- Check if critical Richardson number is exceeded. Further, check
1644	!-- if there is a maxium in the wind profile in order to detect also
1645	!-- ABL heights in the stable boundary layer.
1646	IF ( zi_local == 0.0_wp .AND. &
1647	( ri_bulk > ri_bulk_crit .OR. k == k_max_loc ) ) &
1648	zi_local = zu(k)
1649	ENDDO
1650	!
1651	!-- Assure that the minimum local boundary-layer depth is at least at
1652	!-- the second vertical grid level.
1653	zi_l = zi_l + MAX( zi_local, zu(k_surface+2) )
1654
1655	ENDDO
1656
1657	ENDIF
1658	!
1659	!-- Do the same at the north and south boundaries.
1660	IF ( bc_dirichlet_s .OR. bc_dirichlet_n ) THEN
1661
1662	num_boundary_gp_non_cyclic_l = num_boundary_gp_non_cyclic_l + &
1663	nxr - nxl + 1
1664
1665	j = MERGE( -1, nyn + 1, bc_dirichlet_s )
1666
1667	DO i = nxl, nxr
1668	k_surface = topo_top_ind(j,i,0)
1669
1670	IF ( humidity ) THEN
1671	vpt_surface = pt(k_surface+2,j,i) * &
1672	( 1.0_wp + 0.61_wp * q(k_surface+2,j,i) )
1673	vpt_col = pt(:,j,i) * ( 1.0_wp + 0.61_wp * q(:,j,i) )
1674	ELSE
1675	vpt_surface = pt(k_surface+2,j,i)
1676	vpt_col = pt(:,j,i)
1677	ENDIF
1678
1679	uv_abs(:) = SQRT( u(:,j,i)**2 + &
1680	MERGE( v(:,j+1,i), v(:,j,i), &
1681	bc_dirichlet_s )**2 )
1682	!
1683	!-- Determine index of the maximum wind speed
1684	k_max_loc = MAXLOC( uv_abs(:), DIM = 1 ) - 1
1685
1686	zi_local = 0.0_wp
1687	DO k = k_surface+1, nzt
1688	ri_bulk = zu(k) * g / vpt_surface * &
1689	( vpt_col(k) - vpt_surface ) / &
1690	( uv_abs(k) + 1E-5_wp )
1691	!
1692	!-- Check if critical Richardson number is exceeded. Further, check
1693	!-- if there is a maxium in the wind profile in order to detect also
1694	!-- ABL heights in the stable boundary layer.
1695	IF ( zi_local == 0.0_wp .AND. &
1696	( ri_bulk > ri_bulk_crit .OR. k == k_max_loc ) ) &
1697	zi_local = zu(k)
1698	ENDDO
1699	zi_l = zi_l + MAX( zi_local, zu(k_surface+2) )
1700
1701	ENDDO
1702
1703	ENDIF
1704
1705	#if defined( __parallel )
1706	CALL MPI_ALLREDUCE( zi_l, zi_ribulk, 1, MPI_REAL, MPI_SUM, &
1707	comm2d, ierr )
1708	CALL MPI_ALLREDUCE( num_boundary_gp_non_cyclic_l, &
1709	num_boundary_gp_non_cyclic, &
1710	1, MPI_INTEGER, MPI_SUM, comm2d, ierr )
1711	#else
1712	zi_ribulk = zi_l
1713	num_boundary_gp_non_cyclic = num_boundary_gp_non_cyclic_l
1714	#endif
1715	zi_ribulk = zi_ribulk / REAL( num_boundary_gp_non_cyclic, KIND = wp )
1716	!
1717	!-- Finally, check if boundary layer depth is not below the any topography.
1718	!-- zi_ribulk will be used to adjust rayleigh damping height, i.e. the
1719	!-- lower level of the sponge layer, as well as to adjust the synthetic
1720	!-- turbulence generator accordingly. If Rayleigh damping would be applied
1721	!-- near buildings, etc., this would spoil the simulation results.
1722	topo_max_l = zw(MAXVAL( topo_top_ind(nys:nyn,nxl:nxr,0) ))
1723
1724	#if defined( __parallel )
1725	CALL MPI_ALLREDUCE( topo_max_l, topo_max, 1, MPI_REAL, MPI_MAX, &
1726	comm2d, ierr )
1727	#else
1728	topo_max = topo_max_l
1729	#endif
1730	! zi_ribulk = MAX( zi_ribulk, topo_max )
1731
1732	END SUBROUTINE nesting_offl_calc_zi
1733
1734
1735	!------------------------------------------------------------------------------!
1736	! Description:
1737	!------------------------------------------------------------------------------!
1738	!> Adjust the height where the rayleigh damping starts, i.e. the lower level
1739	!> of the sponge layer.
1740	!------------------------------------------------------------------------------!
1741	SUBROUTINE adjust_sponge_layer
1742
1743	INTEGER(iwp) :: k !< loop index in z-direction
1744
1745	REAL(wp) :: rdh !< updated Rayleigh damping height
1746
1747
1748	IF ( rayleigh_damping_height > 0.0_wp .AND. &
1749	rayleigh_damping_factor > 0.0_wp ) THEN
1750	!
1751	!-- Update Rayleigh-damping height and re-calculate height-depending
1752	!-- damping coefficients.
1753	!-- Assure that rayleigh damping starts well above the boundary layer.
1754	rdh = MIN( MAX( zi_ribulk * 1.3_wp, 10.0_wp * dz(1) ), &
1755	0.8_wp * zu(nzt), rayleigh_damping_height )
1756	!
1757	!-- Update Rayleigh damping factor
1758	DO k = nzb+1, nzt
1759	IF ( zu(k) >= rdh ) THEN
1760	rdf(k) = rayleigh_damping_factor * &
1761	( SIN( pi * 0.5_wp * ( zu(k) - rdh ) &
1762	/ ( zu(nzt) - rdh ) ) &
1763	)**2
1764	ENDIF
1765	ENDDO
1766	rdf_sc = rdf
1767
1768	ENDIF
1769
1770	END SUBROUTINE adjust_sponge_layer
1771
1772	!------------------------------------------------------------------------------!
1773	! Description:
1774	! ------------
1775	!> Performs consistency checks
1776	!------------------------------------------------------------------------------!
1777	SUBROUTINE nesting_offl_check_parameters
1778	!
1779	!-- Check if offline nesting is applied in nested child domain.
1780	IF ( nesting_offline .AND. child_domain ) THEN
1781	message_string = 'Offline nesting is only applicable in root model.'
1782	CALL message( 'offline_nesting_check_parameters', 'PA0622', 1, 2, 0, 6, 0 )
1783	ENDIF
1784
1785	END SUBROUTINE nesting_offl_check_parameters
1786
1787	!------------------------------------------------------------------------------!
1788	! Description:
1789	! ------------
1790	!> Reads the parameter list nesting_offl_parameters
1791	!------------------------------------------------------------------------------!
1792	SUBROUTINE nesting_offl_parin
1793
1794	CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file
1795
1796
1797	NAMELIST /nesting_offl_parameters/ nesting_offline
1798
1799	line = ' '
1800
1801	!
1802	!-- Try to find stg package
1803	REWIND ( 11 )
1804	line = ' '
1805	DO WHILE ( INDEX( line, '&nesting_offl_parameters' ) == 0 )
1806	READ ( 11, '(A)', END=20 ) line
1807	ENDDO
1808	BACKSPACE ( 11 )
1809
1810	!
1811	!-- Read namelist
1812	READ ( 11, nesting_offl_parameters, ERR = 10, END = 20 )
1813
1814	GOTO 20
1815
1816	10 BACKSPACE( 11 )
1817	READ( 11 , '(A)') line
1818	CALL parin_fail_message( 'nesting_offl_parameters', line )
1819
1820	20 CONTINUE
1821
1822
1823	END SUBROUTINE nesting_offl_parin
1824
1825	!------------------------------------------------------------------------------!
1826	! Description:
1827	! ------------
1828	!> Writes information about offline nesting into HEADER file
1829	!------------------------------------------------------------------------------!
1830	SUBROUTINE nesting_offl_header ( io )
1831
1832	INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file
1833
1834	WRITE ( io, 1 )
1835	IF ( nesting_offline ) THEN
1836	WRITE ( io, 3 )
1837	ELSE
1838	WRITE ( io, 2 )
1839	ENDIF
1840
1841	1 FORMAT (//' Offline nesting in COSMO model:'/ &
1842	' -------------------------------'/)
1843	2 FORMAT (' --> No offlince nesting is used (default) ')
1844	3 FORMAT (' --> Offlince nesting is used. Boundary data is read from dynamic input file ')
1845
1846	END SUBROUTINE nesting_offl_header
1847
1848	!------------------------------------------------------------------------------!
1849	! Description:
1850	! ------------
1851	!> Allocate arrays used to read boundary data from NetCDF file and initialize
1852	!> boundary data.
1853	!------------------------------------------------------------------------------!
1854	SUBROUTINE nesting_offl_init
1855
1856	INTEGER(iwp) :: n !< running index for chemical species
1857
1858	!
1859	!-- Get time_utc_init from origin_date_time
1860	CALL get_date_time( 0.0_wp, second_of_day = time_utc_init )
1861
1862	!-- Allocate arrays for geostrophic wind components. Arrays will
1863	!-- incorporate 2 time levels in order to interpolate in between.
1864	ALLOCATE( nest_offl%ug(0:1,1:nzt) )
1865	ALLOCATE( nest_offl%vg(0:1,1:nzt) )
1866	!
1867	!-- Allocate arrays for reading left/right boundary values. Arrays will
1868	!-- incorporate 2 time levels in order to interpolate in between. If the core has
1869	!-- no boundary, allocate a dummy array, in order to enable netcdf parallel
1870	!-- access. Dummy arrays will be allocated with dimension length zero.
1871	IF ( bc_dirichlet_l ) THEN
1872	ALLOCATE( nest_offl%u_left(0:1,nzb+1:nzt,nys:nyn) )
1873	ALLOCATE( nest_offl%v_left(0:1,nzb+1:nzt,nysv:nyn) )
1874	ALLOCATE( nest_offl%w_left(0:1,nzb+1:nzt-1,nys:nyn) )
1875	IF ( humidity ) ALLOCATE( nest_offl%q_left(0:1,nzb+1:nzt,nys:nyn) )
1876	IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_left(0:1,nzb+1:nzt,nys:nyn) )
1877	IF ( air_chemistry ) ALLOCATE( nest_offl%chem_left(0:1,nzb+1:nzt,nys:nyn,&
1878	1:UBOUND( chem_species, 1 )) )
1879	ELSE
1880	ALLOCATE( nest_offl%u_left(1:1,1:1,1:1) )
1881	ALLOCATE( nest_offl%v_left(1:1,1:1,1:1) )
1882	ALLOCATE( nest_offl%w_left(1:1,1:1,1:1) )
1883	IF ( humidity ) ALLOCATE( nest_offl%q_left(1:1,1:1,1:1) )
1884	IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_left(1:1,1:1,1:1) )
1885	IF ( air_chemistry ) ALLOCATE( nest_offl%chem_left(1:1,1:1,1:1, &
1886	1:UBOUND( chem_species, 1 )) )
1887	ENDIF
1888	IF ( bc_dirichlet_r ) THEN
1889	ALLOCATE( nest_offl%u_right(0:1,nzb+1:nzt,nys:nyn) )
1890	ALLOCATE( nest_offl%v_right(0:1,nzb+1:nzt,nysv:nyn) )
1891	ALLOCATE( nest_offl%w_right(0:1,nzb+1:nzt-1,nys:nyn) )
1892	IF ( humidity ) ALLOCATE( nest_offl%q_right(0:1,nzb+1:nzt,nys:nyn) )
1893	IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_right(0:1,nzb+1:nzt,nys:nyn) )
1894	IF ( air_chemistry ) ALLOCATE( nest_offl%chem_right(0:1,nzb+1:nzt,nys:nyn,&
1895	1:UBOUND( chem_species, 1 )) )
1896	ELSE
1897	ALLOCATE( nest_offl%u_right(1:1,1:1,1:1) )
1898	ALLOCATE( nest_offl%v_right(1:1,1:1,1:1) )
1899	ALLOCATE( nest_offl%w_right(1:1,1:1,1:1) )
1900	IF ( humidity ) ALLOCATE( nest_offl%q_right(1:1,1:1,1:1) )
1901	IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_right(1:1,1:1,1:1) )
1902	IF ( air_chemistry ) ALLOCATE( nest_offl%chem_right(1:1,1:1,1:1, &
1903	1:UBOUND( chem_species, 1 )) )
1904	ENDIF
1905	!
1906	!-- Allocate arrays for reading north/south boundary values. Arrays will
1907	!-- incorporate 2 time levels in order to interpolate in between. If the core has
1908	!-- no boundary, allocate a dummy array, in order to enable netcdf parallel
1909	!-- access. Dummy arrays will be allocated with dimension length zero.
1910	IF ( bc_dirichlet_n ) THEN
1911	ALLOCATE( nest_offl%u_north(0:1,nzb+1:nzt,nxlu:nxr) )
1912	ALLOCATE( nest_offl%v_north(0:1,nzb+1:nzt,nxl:nxr) )
1913	ALLOCATE( nest_offl%w_north(0:1,nzb+1:nzt-1,nxl:nxr) )
1914	IF ( humidity ) ALLOCATE( nest_offl%q_north(0:1,nzb+1:nzt,nxl:nxr) )
1915	IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_north(0:1,nzb+1:nzt,nxl:nxr) )
1916	IF ( air_chemistry ) ALLOCATE( nest_offl%chem_north(0:1,nzb+1:nzt,nxl:nxr,&
1917	1:UBOUND( chem_species, 1 )) )
1918	ELSE
1919	ALLOCATE( nest_offl%u_north(1:1,1:1,1:1) )
1920	ALLOCATE( nest_offl%v_north(1:1,1:1,1:1) )
1921	ALLOCATE( nest_offl%w_north(1:1,1:1,1:1) )
1922	IF ( humidity ) ALLOCATE( nest_offl%q_north(1:1,1:1,1:1) )
1923	IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_north(1:1,1:1,1:1) )
1924	IF ( air_chemistry ) ALLOCATE( nest_offl%chem_north(1:1,1:1,1:1, &
1925	1:UBOUND( chem_species, 1 )) )
1926	ENDIF
1927	IF ( bc_dirichlet_s ) THEN
1928	ALLOCATE( nest_offl%u_south(0:1,nzb+1:nzt,nxlu:nxr) )
1929	ALLOCATE( nest_offl%v_south(0:1,nzb+1:nzt,nxl:nxr) )
1930	ALLOCATE( nest_offl%w_south(0:1,nzb+1:nzt-1,nxl:nxr) )
1931	IF ( humidity ) ALLOCATE( nest_offl%q_south(0:1,nzb+1:nzt,nxl:nxr) )
1932	IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_south(0:1,nzb+1:nzt,nxl:nxr) )
1933	IF ( air_chemistry ) ALLOCATE( nest_offl%chem_south(0:1,nzb+1:nzt,nxl:nxr,&
1934	1:UBOUND( chem_species, 1 )) )
1935	ELSE
1936	ALLOCATE( nest_offl%u_south(1:1,1:1,1:1) )
1937	ALLOCATE( nest_offl%v_south(1:1,1:1,1:1) )
1938	ALLOCATE( nest_offl%w_south(1:1,1:1,1:1) )
1939	IF ( humidity ) ALLOCATE( nest_offl%q_south(1:1,1:1,1:1) )
1940	IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_south(1:1,1:1,1:1) )
1941	IF ( air_chemistry ) ALLOCATE( nest_offl%chem_south(1:1,1:1,1:1, &
1942	1:UBOUND( chem_species, 1 )) )
1943	ENDIF
1944	!
1945	!-- Allocate arrays for reading data at the top boundary. In contrast to the
1946	!-- lateral boundaries, every core reads these data so that no dummy
1947	!-- arrays need to be allocated.
1948	ALLOCATE( nest_offl%u_top(0:1,nys:nyn,nxlu:nxr) )
1949	ALLOCATE( nest_offl%v_top(0:1,nysv:nyn,nxl:nxr) )
1950	ALLOCATE( nest_offl%w_top(0:1,nys:nyn,nxl:nxr) )
1951	IF ( humidity ) ALLOCATE( nest_offl%q_top(0:1,nys:nyn,nxl:nxr) )
1952	IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_top(0:1,nys:nyn,nxl:nxr) )
1953	IF ( air_chemistry ) ALLOCATE( nest_offl%chem_top(0:1,nys:nyn,nxl:nxr, &
1954	1:UBOUND( chem_species, 1 )) )
1955	!
1956	!-- For chemical species, create the names of the variables. This is necessary
1957	!-- to identify the respective variable and write it onto the correct array
1958	!-- in the chem_species datatype.
1959	IF ( air_chemistry ) THEN
1960	ALLOCATE( nest_offl%chem_from_file_l(1:UBOUND( chem_species, 1 )) )
1961	ALLOCATE( nest_offl%chem_from_file_n(1:UBOUND( chem_species, 1 )) )
1962	ALLOCATE( nest_offl%chem_from_file_r(1:UBOUND( chem_species, 1 )) )
1963	ALLOCATE( nest_offl%chem_from_file_s(1:UBOUND( chem_species, 1 )) )
1964	ALLOCATE( nest_offl%chem_from_file_t(1:UBOUND( chem_species, 1 )) )
1965
1966	ALLOCATE( nest_offl%var_names_chem_l(1:UBOUND( chem_species, 1 )) )
1967	ALLOCATE( nest_offl%var_names_chem_n(1:UBOUND( chem_species, 1 )) )
1968	ALLOCATE( nest_offl%var_names_chem_r(1:UBOUND( chem_species, 1 )) )
1969	ALLOCATE( nest_offl%var_names_chem_s(1:UBOUND( chem_species, 1 )) )
1970	ALLOCATE( nest_offl%var_names_chem_t(1:UBOUND( chem_species, 1 )) )
1971	!
1972	!-- Initialize flags that indicate whether the variable is on file or
1973	!-- not. Please note, this is only necessary for chemistry variables.
1974	nest_offl%chem_from_file_l(:) = .FALSE.
1975	nest_offl%chem_from_file_n(:) = .FALSE.
1976	nest_offl%chem_from_file_r(:) = .FALSE.
1977	nest_offl%chem_from_file_s(:) = .FALSE.
1978	nest_offl%chem_from_file_t(:) = .FALSE.
1979
1980	DO n = 1, UBOUND( chem_species, 1 )
1981	nest_offl%var_names_chem_l(n) = nest_offl%char_l // &
1982	TRIM(chem_species(n)%name)
1983	nest_offl%var_names_chem_n(n) = nest_offl%char_n // &
1984	TRIM(chem_species(n)%name)
1985	nest_offl%var_names_chem_r(n) = nest_offl%char_r // &
1986	TRIM(chem_species(n)%name)
1987	nest_offl%var_names_chem_s(n) = nest_offl%char_s // &
1988	TRIM(chem_species(n)%name)
1989	nest_offl%var_names_chem_t(n) = nest_offl%char_t // &
1990	TRIM(chem_species(n)%name)
1991	ENDDO
1992	ENDIF
1993	!
1994	!-- Before initial data input is initiated, check if dynamic input file is
1995	!-- present.
1996	IF ( .NOT. input_pids_dynamic ) THEN
1997	message_string = 'nesting_offline = .TRUE. requires dynamic ' // &
1998	'input file ' // &
1999	TRIM( input_file_dynamic ) // TRIM( coupling_char )
2000	CALL message( 'nesting_offl_init', 'PA0546', 1, 2, 0, 6, 0 )
2001	ENDIF
2002	!
2003	!-- Read COSMO data at lateral and top boundaries
2004	CALL nesting_offl_input
2005	!
2006	!-- Check if sufficient time steps are provided to cover the entire
2007	!-- simulation. Note, dynamic input is only required for the 3D simulation,
2008	!-- not for the soil/wall spinup. However, as the spinup time is added
2009	!-- to the end_time, this must be considered here.
2010	IF ( end_time - spinup_time > &
2011	nest_offl%time(nest_offl%nt-1) - time_utc_init ) THEN
2012	message_string = 'end_time of the simulation exceeds the ' // &
2013	'time dimension in the dynamic input file.'
2014	CALL message( 'nesting_offl_init', 'PA0183', 1, 2, 0, 6, 0 )
2015	ENDIF
2016
2017	IF ( nest_offl%time(0) /= time_utc_init ) THEN
2018	message_string = 'Offline nesting: time dimension must start at ' // &
2019	' time_utc_init.'
2020	CALL message( 'nesting_offl_init', 'PA0676', 1, 2, 0, 6, 0 )
2021	ENDIF
2022	!
2023	!-- Initialize boundary data. Please note, do not initialize boundaries in
2024	!-- case of restart runs. This case the boundaries are already initialized
2025	!-- and the boundary data from file would be on the wrong time level.
2026	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
2027	IF ( bc_dirichlet_l ) THEN
2028	u(nzb+1:nzt,nys:nyn,0) = nest_offl%u_left(0,nzb+1:nzt,nys:nyn)
2029	v(nzb+1:nzt,nysv:nyn,-1) = nest_offl%v_left(0,nzb+1:nzt,nysv:nyn)
2030	w(nzb+1:nzt-1,nys:nyn,-1) = nest_offl%w_left(0,nzb+1:nzt-1,nys:nyn)
2031	IF ( .NOT. neutral ) pt(nzb+1:nzt,nys:nyn,-1) = &
2032	nest_offl%pt_left(0,nzb+1:nzt,nys:nyn)
2033	IF ( humidity ) q(nzb+1:nzt,nys:nyn,-1) = &
2034	nest_offl%q_left(0,nzb+1:nzt,nys:nyn)
2035	IF ( air_chemistry ) THEN
2036	DO n = 1, UBOUND( chem_species, 1 )
2037	IF( nest_offl%chem_from_file_l(n) ) THEN
2038	chem_species(n)%conc(nzb+1:nzt,nys:nyn,-1) = &
2039	nest_offl%chem_left(0,nzb+1:nzt,nys:nyn,n)
2040	ENDIF
2041	ENDDO
2042	ENDIF
2043	ENDIF
2044	IF ( bc_dirichlet_r ) THEN
2045	u(nzb+1:nzt,nys:nyn,nxr+1) = nest_offl%u_right(0,nzb+1:nzt,nys:nyn)
2046	v(nzb+1:nzt,nysv:nyn,nxr+1) = nest_offl%v_right(0,nzb+1:nzt,nysv:nyn)
2047	w(nzb+1:nzt-1,nys:nyn,nxr+1) = nest_offl%w_right(0,nzb+1:nzt-1,nys:nyn)
2048	IF ( .NOT. neutral ) pt(nzb+1:nzt,nys:nyn,nxr+1) = &
2049	nest_offl%pt_right(0,nzb+1:nzt,nys:nyn)
2050	IF ( humidity ) q(nzb+1:nzt,nys:nyn,nxr+1) = &
2051	nest_offl%q_right(0,nzb+1:nzt,nys:nyn)
2052	IF ( air_chemistry ) THEN
2053	DO n = 1, UBOUND( chem_species, 1 )
2054	IF( nest_offl%chem_from_file_r(n) ) THEN
2055	chem_species(n)%conc(nzb+1:nzt,nys:nyn,nxr+1) = &
2056	nest_offl%chem_right(0,nzb+1:nzt,nys:nyn,n)
2057	ENDIF
2058	ENDDO
2059	ENDIF
2060	ENDIF
2061	IF ( bc_dirichlet_s ) THEN
2062	u(nzb+1:nzt,-1,nxlu:nxr) = nest_offl%u_south(0,nzb+1:nzt,nxlu:nxr)
2063	v(nzb+1:nzt,0,nxl:nxr) = nest_offl%v_south(0,nzb+1:nzt,nxl:nxr)
2064	w(nzb+1:nzt-1,-1,nxl:nxr) = nest_offl%w_south(0,nzb+1:nzt-1,nxl:nxr)
2065	IF ( .NOT. neutral ) pt(nzb+1:nzt,-1,nxl:nxr) = &
2066	nest_offl%pt_south(0,nzb+1:nzt,nxl:nxr)
2067	IF ( humidity ) q(nzb+1:nzt,-1,nxl:nxr) = &
2068	nest_offl%q_south(0,nzb+1:nzt,nxl:nxr)
2069	IF ( air_chemistry ) THEN
2070	DO n = 1, UBOUND( chem_species, 1 )
2071	IF( nest_offl%chem_from_file_s(n) ) THEN
2072	chem_species(n)%conc(nzb+1:nzt,-1,nxl:nxr) = &
2073	nest_offl%chem_south(0,nzb+1:nzt,nxl:nxr,n)
2074	ENDIF
2075	ENDDO
2076	ENDIF
2077	ENDIF
2078	IF ( bc_dirichlet_n ) THEN
2079	u(nzb+1:nzt,nyn+1,nxlu:nxr) = nest_offl%u_north(0,nzb+1:nzt,nxlu:nxr)
2080	v(nzb+1:nzt,nyn+1,nxl:nxr) = nest_offl%v_north(0,nzb+1:nzt,nxl:nxr)
2081	w(nzb+1:nzt-1,nyn+1,nxl:nxr) = nest_offl%w_north(0,nzb+1:nzt-1,nxl:nxr)
2082	IF ( .NOT. neutral ) pt(nzb+1:nzt,nyn+1,nxl:nxr) = &
2083	nest_offl%pt_north(0,nzb+1:nzt,nxl:nxr)
2084	IF ( humidity ) q(nzb+1:nzt,nyn+1,nxl:nxr) = &
2085	nest_offl%q_north(0,nzb+1:nzt,nxl:nxr)
2086	IF ( air_chemistry ) THEN
2087	DO n = 1, UBOUND( chem_species, 1 )
2088	IF( nest_offl%chem_from_file_n(n) ) THEN
2089	chem_species(n)%conc(nzb+1:nzt,nyn+1,nxl:nxr) = &
2090	nest_offl%chem_north(0,nzb+1:nzt,nxl:nxr,n)
2091	ENDIF
2092	ENDDO
2093	ENDIF
2094	ENDIF
2095	!
2096	!-- Initialize geostrophic wind components. Actually this is already done in
2097	!-- init_3d_model when initializing_action = 'inifor', however, in speical
2098	!-- case of user-defined initialization this will be done here again, in
2099	!-- order to have a consistent initialization.
2100	ug(nzb+1:nzt) = nest_offl%ug(0,nzb+1:nzt)
2101	vg(nzb+1:nzt) = nest_offl%vg(0,nzb+1:nzt)
2102	!
2103	!-- Set bottom and top boundary condition for geostrophic wind components
2104	ug(nzt+1) = ug(nzt)
2105	vg(nzt+1) = vg(nzt)
2106	ug(nzb) = ug(nzb+1)
2107	vg(nzb) = vg(nzb+1)
2108	ENDIF
2109	!
2110	!-- After boundary data is initialized, mask topography at the
2111	!-- boundaries for the velocity components.
2112	u = MERGE( u, 0.0_wp, BTEST( wall_flags_0, 1 ) )
2113	v = MERGE( v, 0.0_wp, BTEST( wall_flags_0, 2 ) )
2114	w = MERGE( w, 0.0_wp, BTEST( wall_flags_0, 3 ) )
2115	!
2116	!-- Initial calculation of the boundary layer depth from the prescribed
2117	!-- boundary data. This is requiered for initialize the synthetic turbulence
2118	!-- generator correctly.
2119	CALL nesting_offl_calc_zi
2120
2121	!
2122	!-- After boundary data is initialized, ensure mass conservation. Not
2123	!-- necessary in restart runs.
2124	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
2125	CALL nesting_offl_mass_conservation
2126	ENDIF
2127
2128	END SUBROUTINE nesting_offl_init
2129
2130	!------------------------------------------------------------------------------!
2131	! Description:
2132	!------------------------------------------------------------------------------!
2133	!> Interpolation function, used to interpolate boundary data in time.
2134	!------------------------------------------------------------------------------!
2135	FUNCTION interpolate_in_time( var_t1, var_t2, fac )
2136
2137	REAL(wp) :: interpolate_in_time !< time-interpolated boundary value
2138	REAL(wp) :: var_t1 !< boundary value at t1
2139	REAL(wp) :: var_t2 !< boundary value at t2
2140	REAL(wp) :: fac !< interpolation factor
2141
2142	interpolate_in_time = ( 1.0_wp - fac ) * var_t1 + fac * var_t2
2143
2144	END FUNCTION interpolate_in_time
2145
2146
2147
2148	END MODULE nesting_offl_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |