Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |