Home

Context Navigation

nesting_offl_mod.f90 @ 4271

Last change on this file since 4271 was 4270, checked in by monakurppa, 5 years ago

Updates for the aerosol module salsa

Reformat the emission data for the pre-processed mode (lod=2) to limit the size of the input file. Now, the emission is given as total number emission per emission category (ncat) and the number size distribution is given in variable emission_number_fracs(ncat,nbins).
Implement offline nesting for salsa. Additionally, offline nesting can be turn off specifically for salsa using the switch nesting_offline_salsa.
Remove init_aerosol_type and init_gases_type from salsa_parin and define them based on the initializing_actions
parameter definition removed from "season" and "season_z01" is added to parin
bugfix in application of index_hh after implementing the new palm_date_time_mod

Property svn:keywords set to Id

File size: 111.3 KB

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

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

Context Navigation

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

Download in other formats: