Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |