Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |