Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |