Home

Context Navigation

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

Last change on this file since 4416 was 4360, checked in by suehring, 5 years ago
Bugfix in output of time-averaged plant-canopy quanities; Output of plant-canopy data only where tall canopy is defined; land-surface model: fix wrong location strings; tests: update urban test case; all source code files: copyright update
Property svn:keywords set to `Id`
File size: 110.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |