Home

Context Navigation

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

Last change on this file since 4282 was 4273, checked in by monakurppa, 5 years ago
Add logical switched nesting_chem and nesting_offline_chem
Property svn:keywords set to `Id`
File size: 110.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |