Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |