Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |