Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |