Home

Context Navigation

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

Last change on this file since 4631 was 4582, checked in by suehring, 4 years ago
remove unused variable from last commit
Property svn:keywords set to `Id`
File size: 111.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |