Home

Context Navigation

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

Last change on this file since 4239 was 4231, checked in by suehring, 5 years ago
Bugfix in array deallocation
Property svn:keywords set to `Id`
File size: 110.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |