Home

Context Navigation

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

Last change on this file since 4561 was 4561, checked in by suehring, 4 years ago
mesoscale nesting: Adapt mass-flux correction also for the anelastic approximation
Property svn:keywords set to `Id`
File size: 111.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |