Home

Context Navigation

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

Last change on this file since 4180 was 4180, checked in by scharf, 5 years ago
removed comments in 'Former revisions' section that are older than 01.01.2019
Property svn:keywords set to `Id`
File size: 65.8 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	!
[3705]	23	!
[3347]	24	! Former revisions:
	25	! -----------------
[3413]	26	! $Id: nesting_offl_mod.f90 4180 2019-08-21 14:37:54Z scharf $
[4169]	27	! Additional check added.
	28	!
	29	! 4168 2019-08-16 13:50:17Z suehring
[4168]	30	! Replace function get_topography_top_index by topo_top_ind
	31	!
	32	! 4125 2019-07-29 13:31:44Z suehring
[4125]	33	! In order to enable netcdf parallel access, allocate dummy arrays for the
	34	! lateral boundary data on cores that actually do not belong to these
	35	! boundaries.
	36	!
	37	! 4079 2019-07-09 18:04:41Z suehring
[4079]	38	! - Set boundary condition for w at nzt+1 at the lateral boundaries, even
	39	! though these won't enter the numerical solution. However, due to the mass
	40	! conservation these values might some up to very large values which will
	41	! occur in the run-control file
	42	! - Bugfix in offline nesting of chemical species
	43	! - Do not set Neumann conditions for TKE and passive scalar
	44	!
	45	! 4022 2019-06-12 11:52:39Z suehring
[4022]	46	! Detection of boundary-layer depth in stable boundary layer on basis of
	47	! boundary data improved
	48	! Routine for boundary-layer depth calculation renamed and made public
	49	!
	50	! 3987 2019-05-22 09:52:13Z kanani
[3987]	51	! Introduce alternative switch for debug output during timestepping
	52	!
	53	! 3964 2019-05-09 09:48:32Z suehring
[3964]	54	! Ensure that veloctiy term in calculation of bulk Richardson number does not
	55	! become zero
	56	!
	57	! 3937 2019-04-29 15:09:07Z suehring
[3937]	58	! Set boundary conditon on upper-left and upper-south grid point for the u- and
	59	! v-component, respectively.
	60	!
	61	! 3891 2019-04-12 17:52:01Z suehring
[3891]	62	! Bugfix, do not overwrite lateral and top boundary data in case of restart
	63	! runs.
	64	!
	65	! 3885 2019-04-11 11:29:34Z kanani
[3885]	66	! Changes related to global restructuring of location messages and introduction
	67	! of additional debug messages
	68	!
	69	!
[3858]	70	! Do local data exchange for chemistry variables only when boundary data is
	71	! coming from dynamic file
	72	!
	73	! 3737 2019-02-12 16:57:06Z suehring
[3737]	74	! Introduce mesoscale nesting for chemical species
	75	!
	76	! 3705 2019-01-29 19:56:39Z suehring
[3705]	77	! Formatting adjustments
	78	!
	79	! 3704 2019-01-29 19:51:41Z suehring
[3579]	80	! Check implemented for offline nesting in child domain
	81	!
[3413]	82	!
[3347]	83	! Description:
	84	! ------------
	85	!> Offline nesting in larger-scale models. Boundary conditions for the simulation
	86	!> are read from NetCDF file and are prescribed onto the respective arrays.
	87	!> Further, a mass-flux correction is performed to maintain the mass balance.
	88	!--------------------------------------------------------------------------------!
	89	MODULE nesting_offl_mod
	90
	91	USE arrays_3d, &
	92	ONLY: dzw, e, diss, pt, pt_init, q, q_init, s, u, u_init, ug, v, &
[3737]	93	v_init, vg, w, zu, zw
	94
[3876]	95	USE chem_modules, &
[3737]	96	ONLY: chem_species
[3347]	97
	98	USE control_parameters, &
[4169]	99	ONLY: air_chemistry, &
	100	bc_dirichlet_l, &
	101	bc_dirichlet_n, &
	102	bc_dirichlet_r, &
	103	bc_dirichlet_s, &
	104	dt_3d, &
	105	dz, &
	106	constant_diffusion, &
	107	child_domain, &
[3987]	108	debug_output_timestep, &
[4169]	109	end_time, &
[3987]	110	humidity, &
	111	initializing_actions, &
[4169]	112	message_string, &
	113	nesting_offline, &
	114	neutral, &
	115	passive_scalar, &
	116	rans_mode, &
	117	rans_tke_e, &
	118	rayleigh_damping_factor, &
	119	rayleigh_damping_height, &
	120	spinup_time, &
	121	time_since_reference_point, &
	122	volume_flow
[3347]	123
	124	USE cpulog, &
	125	ONLY: cpu_log, log_point
	126
	127	USE grid_variables
	128
	129	USE indices, &
	130	ONLY: nbgp, nx, nxl, nxlg, nxlu, nxr, nxrg, ny, nys, &
[4168]	131	nysv, nysg, nyn, nyng, nzb, nz, nzt, &
	132	topo_top_ind, &
	133	wall_flags_0
[3347]	134
	135	USE kinds
	136
	137	USE netcdf_data_input_mod, &
	138	ONLY: nest_offl
	139
	140	USE pegrid
	141
	142	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
	143	!< bulk Richardson number of 0.25
	144
	145	SAVE
	146	PRIVATE
	147	!
	148	!-- Public subroutines
[4022]	149	PUBLIC nesting_offl_bc, &
	150	nesting_offl_calc_zi, &
	151	nesting_offl_check_parameters, &
	152	nesting_offl_header, &
	153	nesting_offl_init, &
	154	nesting_offl_mass_conservation, &
	155	nesting_offl_parin
[3347]	156	!
	157	!-- Public variables
	158	PUBLIC zi_ribulk
	159
	160	INTERFACE nesting_offl_bc
	161	MODULE PROCEDURE nesting_offl_bc
	162	END INTERFACE nesting_offl_bc
	163
[4022]	164	INTERFACE nesting_offl_calc_zi
	165	MODULE PROCEDURE nesting_offl_calc_zi
	166	END INTERFACE nesting_offl_calc_zi
	167
[3579]	168	INTERFACE nesting_offl_check_parameters
	169	MODULE PROCEDURE nesting_offl_check_parameters
	170	END INTERFACE nesting_offl_check_parameters
	171
[3347]	172	INTERFACE nesting_offl_header
	173	MODULE PROCEDURE nesting_offl_header
	174	END INTERFACE nesting_offl_header
	175
	176	INTERFACE nesting_offl_init
	177	MODULE PROCEDURE nesting_offl_init
	178	END INTERFACE nesting_offl_init
	179
	180	INTERFACE nesting_offl_mass_conservation
	181	MODULE PROCEDURE nesting_offl_mass_conservation
	182	END INTERFACE nesting_offl_mass_conservation
	183
	184	INTERFACE nesting_offl_parin
	185	MODULE PROCEDURE nesting_offl_parin
	186	END INTERFACE nesting_offl_parin
	187
	188	CONTAINS
	189
	190
	191	!------------------------------------------------------------------------------!
	192	! Description:
	193	! ------------
	194	!> In this subroutine a constant mass within the model domain is guaranteed.
	195	!> Larger-scale models may be based on a compressible equation system, which is
	196	!> not consistent with PALMs incompressible equation system. In order to avoid
	197	!> a decrease or increase of mass during the simulation, non-divergent flow
	198	!> through the lateral and top boundaries is compensated by the vertical wind
	199	!> component at the top boundary.
	200	!------------------------------------------------------------------------------!
	201	SUBROUTINE nesting_offl_mass_conservation
	202
	203	IMPLICIT NONE
	204
	205	INTEGER(iwp) :: i !< grid index in x-direction
	206	INTEGER(iwp) :: j !< grid index in y-direction
	207	INTEGER(iwp) :: k !< grid index in z-direction
	208
	209	REAL(wp) :: d_area_t !< inverse of the total area of the horizontal model domain
	210	REAL(wp) :: w_correct !< vertical velocity increment required to compensate non-divergent flow through the boundaries
	211	REAL(wp), DIMENSION(1:3) :: volume_flow_l !< local volume flow
	212
[3987]	213
	214	IF ( debug_output_timestep ) CALL debug_message( 'nesting_offl_mass_conservation', 'start' )
	215
[3347]	216	CALL cpu_log( log_point(58), 'offline nesting', 'start' )
	217
	218	volume_flow = 0.0_wp
	219	volume_flow_l = 0.0_wp
	220
	221	d_area_t = 1.0_wp / ( ( nx + 1 ) * dx * ( ny + 1 ) * dy )
	222
	223	IF ( bc_dirichlet_l ) THEN
	224	i = nxl
	225	DO j = nys, nyn
	226	DO k = nzb+1, nzt
	227	volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k) * dy &
	228	* MERGE( 1.0_wp, 0.0_wp, &
	229	BTEST( wall_flags_0(k,j,i), 1 ) )
	230	ENDDO
	231	ENDDO
	232	ENDIF
	233	IF ( bc_dirichlet_r ) THEN
	234	i = nxr+1
	235	DO j = nys, nyn
	236	DO k = nzb+1, nzt
	237	volume_flow_l(1) = volume_flow_l(1) - u(k,j,i) * dzw(k) * dy &
	238	* MERGE( 1.0_wp, 0.0_wp, &
	239	BTEST( wall_flags_0(k,j,i), 1 ) )
	240	ENDDO
	241	ENDDO
	242	ENDIF
	243	IF ( bc_dirichlet_s ) THEN
	244	j = nys
	245	DO i = nxl, nxr
	246	DO k = nzb+1, nzt
	247	volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k) * dx &
	248	* MERGE( 1.0_wp, 0.0_wp, &
	249	BTEST( wall_flags_0(k,j,i), 2 ) )
	250	ENDDO
	251	ENDDO
	252	ENDIF
	253	IF ( bc_dirichlet_n ) THEN
	254	j = nyn+1
	255	DO i = nxl, nxr
	256	DO k = nzb+1, nzt
	257	volume_flow_l(2) = volume_flow_l(2) - v(k,j,i) * dzw(k) * dx &
	258	* MERGE( 1.0_wp, 0.0_wp, &
	259	BTEST( wall_flags_0(k,j,i), 2 ) )
	260	ENDDO
	261	ENDDO
	262	ENDIF
	263	!
	264	!-- Top boundary
	265	k = nzt
	266	DO i = nxl, nxr
	267	DO j = nys, nyn
	268	volume_flow_l(3) = volume_flow_l(3) - w(k,j,i) * dx * dy
	269	ENDDO
	270	ENDDO
	271
	272	#if defined( __parallel )
	273	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
	274	CALL MPI_ALLREDUCE( volume_flow_l, volume_flow, 3, MPI_REAL, MPI_SUM, &
	275	comm2d, ierr )
	276	#else
	277	volume_flow = volume_flow_l
	278	#endif
	279
	280	w_correct = SUM( volume_flow ) * d_area_t
	281
	282	DO i = nxl, nxr
	283	DO j = nys, nyn
	284	DO k = nzt, nzt + 1
	285	w(k,j,i) = w(k,j,i) + w_correct
	286	ENDDO
	287	ENDDO
	288	ENDDO
	289
	290	CALL cpu_log( log_point(58), 'offline nesting', 'stop' )
	291
[3987]	292	IF ( debug_output_timestep ) CALL debug_message( 'nesting_offl_mass_conservation', 'end' )
	293
[3347]	294	END SUBROUTINE nesting_offl_mass_conservation
	295
	296
	297	!------------------------------------------------------------------------------!
	298	! Description:
	299	! ------------
	300	!> Set the lateral and top boundary conditions in case the PALM domain is
	301	!> nested offline in a mesoscale model. Further, average boundary data and
	302	!> determine mean profiles, further used for correct damping in the sponge
	303	!> layer.
	304	!------------------------------------------------------------------------------!
	305	SUBROUTINE nesting_offl_bc
	306
	307	IMPLICIT NONE
	308
	309	INTEGER(iwp) :: i !< running index x-direction
	310	INTEGER(iwp) :: j !< running index y-direction
	311	INTEGER(iwp) :: k !< running index z-direction
[3737]	312	INTEGER(iwp) :: n !< running index for chemical species
[3347]	313
	314	REAL(wp) :: fac_dt !< interpolation factor
	315
	316	REAL(wp), DIMENSION(nzb:nzt+1) :: pt_ref !< reference profile for potential temperature
	317	REAL(wp), DIMENSION(nzb:nzt+1) :: pt_ref_l !< reference profile for potential temperature on subdomain
	318	REAL(wp), DIMENSION(nzb:nzt+1) :: q_ref !< reference profile for mixing ratio on subdomain
	319	REAL(wp), DIMENSION(nzb:nzt+1) :: q_ref_l !< reference profile for mixing ratio on subdomain
	320	REAL(wp), DIMENSION(nzb:nzt+1) :: u_ref !< reference profile for u-component on subdomain
	321	REAL(wp), DIMENSION(nzb:nzt+1) :: u_ref_l !< reference profile for u-component on subdomain
	322	REAL(wp), DIMENSION(nzb:nzt+1) :: v_ref !< reference profile for v-component on subdomain
	323	REAL(wp), DIMENSION(nzb:nzt+1) :: v_ref_l !< reference profile for v-component on subdomain
	324
[3885]	325
[3987]	326	IF ( debug_output_timestep ) CALL debug_message( 'nesting_offl_bc', 'start' )
	327
[3347]	328	CALL cpu_log( log_point(58), 'offline nesting', 'start' )
	329	!
	330	!-- Set mean profiles, derived from boundary data, to zero
	331	pt_ref = 0.0_wp
	332	q_ref = 0.0_wp
	333	u_ref = 0.0_wp
	334	v_ref = 0.0_wp
	335
	336	pt_ref_l = 0.0_wp
	337	q_ref_l = 0.0_wp
	338	u_ref_l = 0.0_wp
	339	v_ref_l = 0.0_wp
	340	!
	341	!-- Determine interpolation factor and limit it to 1. This is because
	342	!-- t+dt can slightly exceed time(tind_p) before boundary data is updated
	343	!-- again.
	344	fac_dt = ( time_since_reference_point - nest_offl%time(nest_offl%tind) &
	345	+ dt_3d ) / &
	346	( nest_offl%time(nest_offl%tind_p) - nest_offl%time(nest_offl%tind) )
	347	fac_dt = MIN( 1.0_wp, fac_dt )
	348	!
	349	!-- Set boundary conditions of u-, v-, w-component, as well as q, and pt.
	350	!-- Note, boundary values at the left boundary: i=-1 (v,w,pt,q) and
	351	!-- i=0 (u), at the right boundary: i=nxr+1 (all), at the south boundary:
	352	!-- j=-1 (u,w,pt,q) and j=0 (v), at the north boundary: j=nyn+1 (all).
	353	!-- Please note, at the left (for u) and south (for v) boundary, values
	354	!-- for u and v are set also at i/j=-1, since these values are used in
	355	!-- boundary_conditions() to restore prognostic values.
	356	!-- Further, sum up data to calculate mean profiles from boundary data,
	357	!-- used for Rayleigh damping.
	358	IF ( bc_dirichlet_l ) THEN
	359
	360	DO j = nys, nyn
	361	DO k = nzb+1, nzt
	362	u(k,j,0) = interpolate_in_time( nest_offl%u_left(0,k,j), &
	363	nest_offl%u_left(1,k,j), &
	364	fac_dt ) * &
	365	MERGE( 1.0_wp, 0.0_wp, &
	366	BTEST( wall_flags_0(k,j,0), 1 ) )
	367	u(k,j,-1) = u(k,j,0)
	368	ENDDO
	369	u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,j,0)
	370	ENDDO
	371
	372	DO j = nys, nyn
	373	DO k = nzb+1, nzt-1
	374	w(k,j,-1) = interpolate_in_time( nest_offl%w_left(0,k,j), &
	375	nest_offl%w_left(1,k,j), &
	376	fac_dt ) * &
	377	MERGE( 1.0_wp, 0.0_wp, &
	378	BTEST( wall_flags_0(k,j,-1), 3 ) )
	379	ENDDO
[4079]	380	w(nzt,j,-1) = w(nzt-1,j,-1)
[3347]	381	ENDDO
	382
	383	DO j = nysv, nyn
	384	DO k = nzb+1, nzt
	385	v(k,j,-1) = interpolate_in_time( nest_offl%v_left(0,k,j), &
	386	nest_offl%v_left(1,k,j), &
	387	fac_dt ) * &
	388	MERGE( 1.0_wp, 0.0_wp, &
	389	BTEST( wall_flags_0(k,j,-1), 2 ) )
	390	ENDDO
	391	v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,j,-1)
	392	ENDDO
	393
	394	IF ( .NOT. neutral ) THEN
	395	DO j = nys, nyn
	396	DO k = nzb+1, nzt
	397	pt(k,j,-1) = interpolate_in_time( nest_offl%pt_left(0,k,j), &
	398	nest_offl%pt_left(1,k,j), &
	399	fac_dt )
	400
	401	ENDDO
	402	pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,j,-1)
	403	ENDDO
	404	ENDIF
	405
	406	IF ( humidity ) THEN
	407	DO j = nys, nyn
	408	DO k = nzb+1, nzt
	409	q(k,j,-1) = interpolate_in_time( nest_offl%q_left(0,k,j), &
	410	nest_offl%q_left(1,k,j), &
	411	fac_dt )
	412
	413	ENDDO
	414	q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,j,-1)
	415	ENDDO
	416	ENDIF
[3737]	417
	418	IF ( air_chemistry ) THEN
	419	DO n = 1, UBOUND( chem_species, 1 )
	420	IF ( nest_offl%chem_from_file_l(n) ) THEN
	421	DO j = nys, nyn
	422	DO k = nzb+1, nzt
	423	chem_species(n)%conc(k,j,-1) = interpolate_in_time( &
	424	nest_offl%chem_left(0,k,j,n),&
	425	nest_offl%chem_left(1,k,j,n),&
	426	fac_dt )
	427	ENDDO
	428	ENDDO
	429	ENDIF
	430	ENDDO
	431	ENDIF
[3347]	432
	433	ENDIF
	434
	435	IF ( bc_dirichlet_r ) THEN
	436
	437	DO j = nys, nyn
	438	DO k = nzb+1, nzt
	439	u(k,j,nxr+1) = interpolate_in_time( nest_offl%u_right(0,k,j), &
	440	nest_offl%u_right(1,k,j), &
	441	fac_dt ) * &
	442	MERGE( 1.0_wp, 0.0_wp, &
	443	BTEST( wall_flags_0(k,j,nxr+1), 1 ) )
	444	ENDDO
	445	u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,j,nxr+1)
	446	ENDDO
	447	DO j = nys, nyn
	448	DO k = nzb+1, nzt-1
	449	w(k,j,nxr+1) = interpolate_in_time( nest_offl%w_right(0,k,j), &
	450	nest_offl%w_right(1,k,j), &
	451	fac_dt ) * &
	452	MERGE( 1.0_wp, 0.0_wp, &
	453	BTEST( wall_flags_0(k,j,nxr+1), 3 ) )
	454	ENDDO
[4079]	455	w(nzt,j,nxr+1) = w(nzt-1,j,nxr+1)
[3347]	456	ENDDO
	457
	458	DO j = nysv, nyn
	459	DO k = nzb+1, nzt
	460	v(k,j,nxr+1) = interpolate_in_time( nest_offl%v_right(0,k,j), &
	461	nest_offl%v_right(1,k,j), &
	462	fac_dt ) * &
	463	MERGE( 1.0_wp, 0.0_wp, &
	464	BTEST( wall_flags_0(k,j,nxr+1), 2 ) )
	465	ENDDO
	466	v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,j,nxr+1)
	467	ENDDO
	468
	469	IF ( .NOT. neutral ) THEN
	470	DO j = nys, nyn
	471	DO k = nzb+1, nzt
	472	pt(k,j,nxr+1) = interpolate_in_time( &
	473	nest_offl%pt_right(0,k,j), &
	474	nest_offl%pt_right(1,k,j), &
	475	fac_dt )
	476	ENDDO
	477	pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,j,nxr+1)
	478	ENDDO
	479	ENDIF
	480
	481	IF ( humidity ) THEN
	482	DO j = nys, nyn
	483	DO k = nzb+1, nzt
	484	q(k,j,nxr+1) = interpolate_in_time( &
	485	nest_offl%q_right(0,k,j), &
	486	nest_offl%q_right(1,k,j), &
	487	fac_dt )
	488
	489	ENDDO
	490	q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,j,nxr+1)
	491	ENDDO
	492	ENDIF
[3737]	493
	494	IF ( air_chemistry ) THEN
	495	DO n = 1, UBOUND( chem_species, 1 )
	496	IF ( nest_offl%chem_from_file_r(n) ) THEN
	497	DO j = nys, nyn
	498	DO k = nzb+1, nzt
	499	chem_species(n)%conc(k,j,nxr+1) = interpolate_in_time(&
	500	nest_offl%chem_right(0,k,j,n),&
	501	nest_offl%chem_right(1,k,j,n),&
	502	fac_dt )
	503	ENDDO
	504	ENDDO
	505	ENDIF
	506	ENDDO
	507	ENDIF
[3347]	508
	509	ENDIF
	510
	511	IF ( bc_dirichlet_s ) THEN
	512
	513	DO i = nxl, nxr
	514	DO k = nzb+1, nzt
	515	v(k,0,i) = interpolate_in_time( nest_offl%v_south(0,k,i), &
	516	nest_offl%v_south(1,k,i), &
	517	fac_dt ) * &
	518	MERGE( 1.0_wp, 0.0_wp, &
	519	BTEST( wall_flags_0(k,0,i), 2 ) )
	520	v(k,-1,i) = v(k,0,i)
	521	ENDDO
	522	v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,0,i)
	523	ENDDO
	524
	525	DO i = nxl, nxr
	526	DO k = nzb+1, nzt-1
	527	w(k,-1,i) = interpolate_in_time( nest_offl%w_south(0,k,i), &
	528	nest_offl%w_south(1,k,i), &
	529	fac_dt ) * &
	530	MERGE( 1.0_wp, 0.0_wp, &
	531	BTEST( wall_flags_0(k,-1,i), 3 ) )
	532	ENDDO
[4079]	533	w(nzt,-1,i) = w(nzt-1,-1,i)
[3347]	534	ENDDO
	535
	536	DO i = nxlu, nxr
	537	DO k = nzb+1, nzt
	538	u(k,-1,i) = interpolate_in_time( nest_offl%u_south(0,k,i), &
	539	nest_offl%u_south(1,k,i), &
	540	fac_dt ) * &
	541	MERGE( 1.0_wp, 0.0_wp, &
	542	BTEST( wall_flags_0(k,-1,i), 1 ) )
	543	ENDDO
	544	u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,-1,i)
	545	ENDDO
	546
	547	IF ( .NOT. neutral ) THEN
	548	DO i = nxl, nxr
	549	DO k = nzb+1, nzt
	550	pt(k,-1,i) = interpolate_in_time( &
	551	nest_offl%pt_south(0,k,i), &
	552	nest_offl%pt_south(1,k,i), &
	553	fac_dt )
	554
	555	ENDDO
	556	pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,-1,i)
	557	ENDDO
	558	ENDIF
	559
	560	IF ( humidity ) THEN
	561	DO i = nxl, nxr
	562	DO k = nzb+1, nzt
	563	q(k,-1,i) = interpolate_in_time( &
	564	nest_offl%q_south(0,k,i), &
	565	nest_offl%q_south(1,k,i), &
	566	fac_dt )
	567
	568	ENDDO
	569	q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,-1,i)
	570	ENDDO
	571	ENDIF
[3737]	572
	573	IF ( air_chemistry ) THEN
	574	DO n = 1, UBOUND( chem_species, 1 )
	575	IF ( nest_offl%chem_from_file_s(n) ) THEN
	576	DO i = nxl, nxr
	577	DO k = nzb+1, nzt
	578	chem_species(n)%conc(k,-1,i) = interpolate_in_time( &
	579	nest_offl%chem_south(0,k,i,n),&
	580	nest_offl%chem_south(1,k,i,n),&
	581	fac_dt )
	582	ENDDO
	583	ENDDO
	584	ENDIF
	585	ENDDO
	586	ENDIF
[3347]	587
	588	ENDIF
	589
	590	IF ( bc_dirichlet_n ) THEN
	591
	592	DO i = nxl, nxr
	593	DO k = nzb+1, nzt
	594	v(k,nyn+1,i) = interpolate_in_time( nest_offl%v_north(0,k,i), &
	595	nest_offl%v_north(1,k,i), &
	596	fac_dt ) * &
	597	MERGE( 1.0_wp, 0.0_wp, &
	598	BTEST( wall_flags_0(k,nyn+1,i), 2 ) )
	599	ENDDO
	600	v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,nyn+1,i)
	601	ENDDO
	602	DO i = nxl, nxr
	603	DO k = nzb+1, nzt-1
	604	w(k,nyn+1,i) = interpolate_in_time( nest_offl%w_north(0,k,i), &
	605	nest_offl%w_north(1,k,i), &
	606	fac_dt ) * &
	607	MERGE( 1.0_wp, 0.0_wp, &
	608	BTEST( wall_flags_0(k,nyn+1,i), 3 ) )
	609	ENDDO
[4079]	610	w(nzt,nyn+1,i) = w(nzt-1,nyn+1,i)
[3347]	611	ENDDO
	612
	613	DO i = nxlu, nxr
	614	DO k = nzb+1, nzt
	615	u(k,nyn+1,i) = interpolate_in_time( nest_offl%u_north(0,k,i), &
	616	nest_offl%u_north(1,k,i), &
	617	fac_dt ) * &
	618	MERGE( 1.0_wp, 0.0_wp, &
	619	BTEST( wall_flags_0(k,nyn+1,i), 1 ) )
	620
	621	ENDDO
	622	u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,nyn+1,i)
	623	ENDDO
	624
	625	IF ( .NOT. neutral ) THEN
	626	DO i = nxl, nxr
	627	DO k = nzb+1, nzt
	628	pt(k,nyn+1,i) = interpolate_in_time( &
	629	nest_offl%pt_north(0,k,i), &
	630	nest_offl%pt_north(1,k,i), &
	631	fac_dt )
	632
	633	ENDDO
	634	pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,nyn+1,i)
	635	ENDDO
	636	ENDIF
	637
	638	IF ( humidity ) THEN
	639	DO i = nxl, nxr
	640	DO k = nzb+1, nzt
	641	q(k,nyn+1,i) = interpolate_in_time( &
	642	nest_offl%q_north(0,k,i), &
	643	nest_offl%q_north(1,k,i), &
	644	fac_dt )
	645
	646	ENDDO
	647	q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,nyn+1,i)
	648	ENDDO
	649	ENDIF
[3737]	650
	651	IF ( air_chemistry ) THEN
	652	DO n = 1, UBOUND( chem_species, 1 )
	653	IF ( nest_offl%chem_from_file_n(n) ) THEN
	654	DO i = nxl, nxr
	655	DO k = nzb+1, nzt
	656	chem_species(n)%conc(k,nyn+1,i) = interpolate_in_time(&
	657	nest_offl%chem_north(0,k,i,n),&
	658	nest_offl%chem_north(1,k,i,n),&
	659	fac_dt )
	660	ENDDO
	661	ENDDO
	662	ENDIF
	663	ENDDO
	664	ENDIF
[3347]	665
	666	ENDIF
	667	!
	668	!-- Top boundary
	669	DO i = nxlu, nxr
	670	DO j = nys, nyn
	671	u(nzt+1,j,i) = interpolate_in_time( nest_offl%u_top(0,j,i), &
	672	nest_offl%u_top(1,j,i), &
	673	fac_dt ) * &
	674	MERGE( 1.0_wp, 0.0_wp, &
	675	BTEST( wall_flags_0(nzt+1,j,i), 1 ) )
	676	u_ref_l(nzt+1) = u_ref_l(nzt+1) + u(nzt+1,j,i)
	677	ENDDO
	678	ENDDO
[3937]	679	!
	680	!-- For left boundary set boundary condition for u-component also at top
	681	!-- grid point.
	682	!-- Note, this has no effect on the numeric solution, only for data output.
	683	IF ( bc_dirichlet_l ) u(nzt+1,:,nxl) = u(nzt+1,:,nxlu)
[3347]	684
	685	DO i = nxl, nxr
	686	DO j = nysv, nyn
	687	v(nzt+1,j,i) = interpolate_in_time( nest_offl%v_top(0,j,i), &
	688	nest_offl%v_top(1,j,i), &
	689	fac_dt ) * &
	690	MERGE( 1.0_wp, 0.0_wp, &
	691	BTEST( wall_flags_0(nzt+1,j,i), 2 ) )
	692	v_ref_l(nzt+1) = v_ref_l(nzt+1) + v(nzt+1,j,i)
	693	ENDDO
	694	ENDDO
[3937]	695	!
	696	!-- For south boundary set boundary condition for v-component also at top
	697	!-- grid point.
	698	!-- Note, this has no effect on the numeric solution, only for data output.
	699	IF ( bc_dirichlet_s ) v(nzt+1,nys,:) = v(nzt+1,nysv,:)
[3347]	700
	701	DO i = nxl, nxr
	702	DO j = nys, nyn
	703	w(nzt,j,i) = interpolate_in_time( nest_offl%w_top(0,j,i), &
	704	nest_offl%w_top(1,j,i), &
	705	fac_dt ) * &
	706	MERGE( 1.0_wp, 0.0_wp, &
	707	BTEST( wall_flags_0(nzt,j,i), 3 ) )
	708	w(nzt+1,j,i) = w(nzt,j,i)
	709	ENDDO
	710	ENDDO
	711
	712
	713	IF ( .NOT. neutral ) THEN
	714	DO i = nxl, nxr
	715	DO j = nys, nyn
	716	pt(nzt+1,j,i) = interpolate_in_time( nest_offl%pt_top(0,j,i), &
	717	nest_offl%pt_top(1,j,i), &
	718	fac_dt )
	719	pt_ref_l(nzt+1) = pt_ref_l(nzt+1) + pt(nzt+1,j,i)
	720	ENDDO
	721	ENDDO
	722	ENDIF
	723
	724	IF ( humidity ) THEN
	725	DO i = nxl, nxr
	726	DO j = nys, nyn
	727	q(nzt+1,j,i) = interpolate_in_time( nest_offl%q_top(0,j,i), &
	728	nest_offl%q_top(1,j,i), &
	729	fac_dt )
	730	q_ref_l(nzt+1) = q_ref_l(nzt+1) + q(nzt+1,j,i)
	731	ENDDO
	732	ENDDO
	733	ENDIF
[3737]	734
	735	IF ( air_chemistry ) THEN
	736	DO n = 1, UBOUND( chem_species, 1 )
	737	IF ( nest_offl%chem_from_file_t(n) ) THEN
	738	DO i = nxl, nxr
	739	DO j = nys, nyn
	740	chem_species(n)%conc(nzt+1,j,i) = interpolate_in_time( &
[4079]	741	nest_offl%chem_top(0,j,i,n), &
	742	nest_offl%chem_top(1,j,i,n), &
[3737]	743	fac_dt )
	744	ENDDO
	745	ENDDO
	746	ENDIF
	747	ENDDO
	748	ENDIF
[3347]	749	!
	750	!-- Moreover, set Neumann boundary condition for subgrid-scale TKE,
	751	!-- passive scalar, dissipation, and chemical species if required
	752	IF ( rans_mode .AND. rans_tke_e ) THEN
	753	IF ( bc_dirichlet_l ) diss(:,:,nxl-1) = diss(:,:,nxl)
	754	IF ( bc_dirichlet_r ) diss(:,:,nxr+1) = diss(:,:,nxr)
	755	IF ( bc_dirichlet_s ) diss(:,nys-1,:) = diss(:,nys,:)
	756	IF ( bc_dirichlet_n ) diss(:,nyn+1,:) = diss(:,nyn,:)
	757	ENDIF
[4079]	758	! IF ( .NOT. constant_diffusion ) THEN
	759	! IF ( bc_dirichlet_l ) e(:,:,nxl-1) = e(:,:,nxl)
	760	! IF ( bc_dirichlet_r ) e(:,:,nxr+1) = e(:,:,nxr)
	761	! IF ( bc_dirichlet_s ) e(:,nys-1,:) = e(:,nys,:)
	762	! IF ( bc_dirichlet_n ) e(:,nyn+1,:) = e(:,nyn,:)
	763	! e(nzt+1,:,:) = e(nzt,:,:)
	764	! ENDIF
	765	! IF ( passive_scalar ) THEN
	766	! IF ( bc_dirichlet_l ) s(:,:,nxl-1) = s(:,:,nxl)
	767	! IF ( bc_dirichlet_r ) s(:,:,nxr+1) = s(:,:,nxr)
	768	! IF ( bc_dirichlet_s ) s(:,nys-1,:) = s(:,nys,:)
	769	! IF ( bc_dirichlet_n ) s(:,nyn+1,:) = s(:,nyn,:)
	770	! ENDIF
[3347]	771
	772	CALL exchange_horiz( u, nbgp )
	773	CALL exchange_horiz( v, nbgp )
	774	CALL exchange_horiz( w, nbgp )
	775	IF ( .NOT. neutral ) CALL exchange_horiz( pt, nbgp )
	776	IF ( humidity ) CALL exchange_horiz( q, nbgp )
[3737]	777	IF ( air_chemistry ) THEN
	778	DO n = 1, UBOUND( chem_species, 1 )
[3858]	779	!
	780	!-- Do local exchange only when necessary, i.e. when data is coming
	781	!-- from dynamic file.
	782	IF ( nest_offl%chem_from_file_t(n) ) &
	783	CALL exchange_horiz( chem_species(n)%conc, nbgp )
[3737]	784	ENDDO
	785	ENDIF
[3347]	786	!
[4079]	787	!-- Set top boundary condition at all horizontal grid points, also at the
	788	!-- lateral boundary grid points.
	789	w(nzt+1,:,:) = w(nzt,:,:)
	790	!
[3347]	791	!-- In case of Rayleigh damping, where the profiles u_init, v_init
	792	!-- q_init and pt_init are still used, update these profiles from the
	793	!-- averaged boundary data.
	794	!-- But first, average these data.
	795	#if defined( __parallel )
	796	CALL MPI_ALLREDUCE( u_ref_l, u_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM, &
	797	comm2d, ierr )
	798	CALL MPI_ALLREDUCE( v_ref_l, v_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM, &
	799	comm2d, ierr )
	800	IF ( humidity ) THEN
	801	CALL MPI_ALLREDUCE( q_ref_l, q_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM, &
	802	comm2d, ierr )
	803	ENDIF
	804	IF ( .NOT. neutral ) THEN
	805	CALL MPI_ALLREDUCE( pt_ref_l, pt_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM,&
	806	comm2d, ierr )
	807	ENDIF
	808	#else
[3704]	809	u_ref = u_ref_l
	810	v_ref = v_ref_l
	811	IF ( humidity ) q_ref = q_ref_l
	812	IF ( .NOT. neutral ) pt_ref = pt_ref_l
[3347]	813	#endif
	814	!
[3704]	815	!-- Average data. Note, reference profiles up to nzt are derived from lateral
	816	!-- boundaries, at the model top it is derived from the top boundary. Thus,
	817	!-- number of input data is different from nzb:nzt compared to nzt+1.
	818	!-- Derived from lateral boundaries.
	819	u_ref(nzb:nzt) = u_ref(nzb:nzt) / REAL( 2.0_wp * ( ny + 1 + nx ), &
	820	KIND = wp )
	821	v_ref(nzb:nzt) = v_ref(nzb:nzt) / REAL( 2.0_wp * ( ny + nx + 1 ), &
	822	KIND = wp )
	823	IF ( humidity ) &
	824	q_ref(nzb:nzt) = q_ref(nzb:nzt) / REAL( 2.0_wp * &
	825	( ny + 1 + nx + 1 ), &
	826	KIND = wp )
	827	IF ( .NOT. neutral ) &
	828	pt_ref(nzb:nzt) = pt_ref(nzb:nzt) / REAL( 2.0_wp * &
	829	( ny + 1 + nx + 1 ), &
	830	KIND = wp )
[3347]	831	!
[3704]	832	!-- Derived from top boundary.
	833	u_ref(nzt+1) = u_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx ), KIND = wp )
	834	v_ref(nzt+1) = v_ref(nzt+1) / REAL( ( ny ) * ( nx + 1 ), KIND = wp )
	835	IF ( humidity ) &
	836	q_ref(nzt+1) = q_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx + 1 ), &
	837	KIND = wp )
	838	IF ( .NOT. neutral ) &
	839	pt_ref(nzt+1) = pt_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx + 1 ), &
	840	KIND = wp )
[3347]	841	!
[3704]	842	!-- Write onto init profiles, which are used for damping
	843	u_init = u_ref
	844	v_init = v_ref
	845	IF ( humidity ) q_init = q_ref
	846	IF ( .NOT. neutral ) pt_init = pt_ref
[3347]	847	!
[3704]	848	!-- Set bottom boundary condition
	849	IF ( humidity ) q_init(nzb) = q_init(nzb+1)
	850	IF ( .NOT. neutral ) pt_init(nzb) = pt_init(nzb+1)
[3347]	851	!
[3704]	852	!-- Further, adjust Rayleigh damping height in case of time-changing conditions.
	853	!-- Therefore, calculate boundary-layer depth first.
[4022]	854	CALL nesting_offl_calc_zi
[3704]	855	CALL adjust_sponge_layer
[3347]	856
	857	!
[3704]	858	!-- Update geostrophic wind components from dynamic input file.
	859	DO k = nzb+1, nzt
	860	ug(k) = interpolate_in_time( nest_offl%ug(0,k), nest_offl%ug(1,k), &
	861	fac_dt )
	862	vg(k) = interpolate_in_time( nest_offl%vg(0,k), nest_offl%vg(1,k), &
	863	fac_dt )
	864	ENDDO
	865	ug(nzt+1) = ug(nzt)
	866	vg(nzt+1) = vg(nzt)
[3347]	867
[3704]	868	CALL cpu_log( log_point(58), 'offline nesting', 'stop' )
[3347]	869
[3987]	870	IF ( debug_output_timestep ) CALL debug_message( 'nesting_offl_bc', 'end' )
	871
	872
[3347]	873	END SUBROUTINE nesting_offl_bc
	874
	875	!------------------------------------------------------------------------------!
	876	! Description:
	877	!------------------------------------------------------------------------------!
	878	!> Calculates the boundary-layer depth from the boundary data, according to
	879	!> bulk-Richardson criterion.
	880	!------------------------------------------------------------------------------!
[4022]	881	SUBROUTINE nesting_offl_calc_zi
[3347]	882
	883	USE basic_constants_and_equations_mod, &
	884	ONLY: g
	885
	886	USE kinds
	887
	888	IMPLICIT NONE
	889
[4022]	890	INTEGER(iwp) :: i !< loop index in x-direction
	891	INTEGER(iwp) :: j !< loop index in y-direction
	892	INTEGER(iwp) :: k !< loop index in z-direction
	893	INTEGER(iwp) :: k_max_loc !< index of maximum wind speed along z-direction
	894	INTEGER(iwp) :: k_surface !< topography top index in z-direction
	895	INTEGER(iwp) :: num_boundary_gp_non_cyclic !< number of non-cyclic boundaries, used for averaging ABL depth
	896	INTEGER(iwp) :: num_boundary_gp_non_cyclic_l !< number of non-cyclic boundaries, used for averaging ABL depth
[3347]	897
	898	REAL(wp) :: ri_bulk !< bulk Richardson number
	899	REAL(wp) :: ri_bulk_crit = 0.25_wp !< critical bulk Richardson number
	900	REAL(wp) :: topo_max !< maximum topography level in model domain
	901	REAL(wp) :: topo_max_l !< maximum topography level in subdomain
	902	REAL(wp) :: vpt_surface !< near-surface virtual potential temperature
	903	REAL(wp) :: zi_l !< mean boundary-layer depth on subdomain
	904	REAL(wp) :: zi_local !< local boundary-layer depth
	905
	906	REAL(wp), DIMENSION(nzb:nzt+1) :: vpt_col !< vertical profile of virtual potential temperature at (j,i)-grid point
[4022]	907	REAL(wp), DIMENSION(nzb:nzt+1) :: uv_abs !< vertical profile of horizontal wind speed at (j,i)-grid point
[3347]	908
	909
	910	!
	911	!-- Calculate mean boundary-layer height from boundary data.
	912	!-- Start with the left and right boundaries.
	913	zi_l = 0.0_wp
[4022]	914	num_boundary_gp_non_cyclic_l = 0
[3347]	915	IF ( bc_dirichlet_l .OR. bc_dirichlet_r ) THEN
	916	!
[4022]	917	!-- Sum-up and store number of boundary grid points used for averaging
	918	!-- ABL depth
	919	num_boundary_gp_non_cyclic_l = num_boundary_gp_non_cyclic_l + &
	920	nxr - nxl + 1
	921	!
[3347]	922	!-- Determine index along x. Please note, index indicates boundary
	923	!-- grid point for scalars.
	924	i = MERGE( -1, nxr + 1, bc_dirichlet_l )
	925
	926	DO j = nys, nyn
	927	!
	928	!-- Determine topography top index at current (j,i) index
[4168]	929	k_surface = topo_top_ind(j,i,0)
[3347]	930	!
	931	!-- Pre-compute surface virtual temperature. Therefore, use 2nd
	932	!-- prognostic level according to Heinze et al. (2017).
	933	IF ( humidity ) THEN
	934	vpt_surface = pt(k_surface+2,j,i) * &
	935	( 1.0_wp + 0.61_wp * q(k_surface+2,j,i) )
	936	vpt_col = pt(:,j,i) * ( 1.0_wp + 0.61_wp * q(:,j,i) )
	937	ELSE
	938	vpt_surface = pt(k_surface+2,j,i)
	939	vpt_col = pt(:,j,i)
	940	ENDIF
	941	!
	942	!-- Calculate local boundary layer height from bulk Richardson number,
	943	!-- i.e. the height where the bulk Richardson number exceeds its
	944	!-- critical value of 0.25 (according to Heinze et al., 2017).
	945	!-- Note, no interpolation of u- and v-component is made, as both
	946	!-- are mainly mean inflow profiles with very small spatial variation.
[3964]	947	!-- Add a safety factor in case the velocity term becomes zero. This
	948	!-- may happen if overhanging 3D structures are directly located at
	949	!-- the boundary, where velocity inside the building is zero
	950	!-- (k_surface is the index of the lowest upward-facing surface).
[4022]	951	uv_abs(:) = SQRT( MERGE( u(:,j,i+1), u(:,j,i), &
	952	bc_dirichlet_l )**2 + &
	953	v(:,j,i)**2 )
	954	!
	955	!-- Determine index of the maximum wind speed
	956	k_max_loc = MAXLOC( uv_abs(:), DIM = 1 ) - 1
	957
[3347]	958	zi_local = 0.0_wp
	959	DO k = k_surface+1, nzt
	960	ri_bulk = zu(k) * g / vpt_surface * &
	961	( vpt_col(k) - vpt_surface ) / &
[4022]	962	( uv_abs(k) + 1E-5_wp )
	963	!
	964	!-- Check if critical Richardson number is exceeded. Further, check
	965	!-- if there is a maxium in the wind profile in order to detect also
	966	!-- ABL heights in the stable boundary layer.
	967	IF ( zi_local == 0.0_wp .AND. &
	968	( ri_bulk > ri_bulk_crit .OR. k == k_max_loc ) ) &
[3347]	969	zi_local = zu(k)
	970	ENDDO
	971	!
	972	!-- Assure that the minimum local boundary-layer depth is at least at
	973	!-- the second vertical grid level.
	974	zi_l = zi_l + MAX( zi_local, zu(k_surface+2) )
	975
	976	ENDDO
	977
	978	ENDIF
	979	!
	980	!-- Do the same at the north and south boundaries.
	981	IF ( bc_dirichlet_s .OR. bc_dirichlet_n ) THEN
	982
[4022]	983	num_boundary_gp_non_cyclic_l = num_boundary_gp_non_cyclic_l + &
	984	nxr - nxl + 1
	985
[3347]	986	j = MERGE( -1, nyn + 1, bc_dirichlet_s )
	987
	988	DO i = nxl, nxr
[4168]	989	k_surface = topo_top_ind(j,i,0)
[3347]	990
	991	IF ( humidity ) THEN
	992	vpt_surface = pt(k_surface+2,j,i) * &
	993	( 1.0_wp + 0.61_wp * q(k_surface+2,j,i) )
	994	vpt_col = pt(:,j,i) * ( 1.0_wp + 0.61_wp * q(:,j,i) )
	995	ELSE
	996	vpt_surface = pt(k_surface+2,j,i)
	997	vpt_col = pt(:,j,i)
	998	ENDIF
	999
[4022]	1000	uv_abs(:) = SQRT( u(:,j,i)**2 + &
	1001	MERGE( v(:,j+1,i), v(:,j,i), &
	1002	bc_dirichlet_s )**2 )
	1003	!
	1004	!-- Determine index of the maximum wind speed
	1005	k_max_loc = MAXLOC( uv_abs(:), DIM = 1 ) - 1
	1006
[3347]	1007	zi_local = 0.0_wp
	1008	DO k = k_surface+1, nzt
	1009	ri_bulk = zu(k) * g / vpt_surface * &
	1010	( vpt_col(k) - vpt_surface ) / &
[4022]	1011	( uv_abs(k) + 1E-5_wp )
	1012	!
	1013	!-- Check if critical Richardson number is exceeded. Further, check
	1014	!-- if there is a maxium in the wind profile in order to detect also
	1015	!-- ABL heights in the stable boundary layer.
	1016	IF ( zi_local == 0.0_wp .AND. &
	1017	( ri_bulk > ri_bulk_crit .OR. k == k_max_loc ) ) &
[3347]	1018	zi_local = zu(k)
	1019	ENDDO
	1020	zi_l = zi_l + MAX( zi_local, zu(k_surface+2) )
	1021
	1022	ENDDO
	1023
	1024	ENDIF
	1025
	1026	#if defined( __parallel )
	1027	CALL MPI_ALLREDUCE( zi_l, zi_ribulk, 1, MPI_REAL, MPI_SUM, &
	1028	comm2d, ierr )
[4022]	1029	CALL MPI_ALLREDUCE( num_boundary_gp_non_cyclic_l, &
	1030	num_boundary_gp_non_cyclic, &
	1031	1, MPI_INTEGER, MPI_SUM, comm2d, ierr )
[3347]	1032	#else
	1033	zi_ribulk = zi_l
[4022]	1034	num_boundary_gp_non_cyclic = num_boundary_gp_non_cyclic_l
[3347]	1035	#endif
[4022]	1036	zi_ribulk = zi_ribulk / REAL( num_boundary_gp_non_cyclic, KIND = wp )
[3347]	1037	!
	1038	!-- Finally, check if boundary layer depth is not below the any topography.
	1039	!-- zi_ribulk will be used to adjust rayleigh damping height, i.e. the
	1040	!-- lower level of the sponge layer, as well as to adjust the synthetic
	1041	!-- turbulence generator accordingly. If Rayleigh damping would be applied
	1042	!-- near buildings, etc., this would spoil the simulation results.
[4168]	1043	topo_max_l = zw(MAXVAL( topo_top_ind(nys:nyn,nxl:nxr,0) ))
[3347]	1044
	1045	#if defined( __parallel )
[4022]	1046	CALL MPI_ALLREDUCE( topo_max_l, topo_max, 1, MPI_REAL, MPI_MAX, &
[3347]	1047	comm2d, ierr )
	1048	#else
	1049	topo_max = topo_max_l
	1050	#endif
[4022]	1051	! zi_ribulk = MAX( zi_ribulk, topo_max )
[3937]	1052
[4022]	1053	END SUBROUTINE nesting_offl_calc_zi
[3347]	1054
	1055
	1056	!------------------------------------------------------------------------------!
	1057	! Description:
	1058	!------------------------------------------------------------------------------!
	1059	!> Adjust the height where the rayleigh damping starts, i.e. the lower level
	1060	!> of the sponge layer.
	1061	!------------------------------------------------------------------------------!
	1062	SUBROUTINE adjust_sponge_layer
	1063
	1064	USE arrays_3d, &
	1065	ONLY: rdf, rdf_sc, zu
	1066
	1067	USE basic_constants_and_equations_mod, &
	1068	ONLY: pi
	1069
	1070	USE kinds
	1071
	1072	IMPLICIT NONE
	1073
	1074	INTEGER(iwp) :: k !< loop index in z-direction
	1075
	1076	REAL(wp) :: rdh !< updated Rayleigh damping height
	1077
	1078
	1079	IF ( rayleigh_damping_height > 0.0_wp .AND. &
	1080	rayleigh_damping_factor > 0.0_wp ) THEN
	1081	!
	1082	!-- Update Rayleigh-damping height and re-calculate height-depending
	1083	!-- damping coefficients.
	1084	!-- Assure that rayleigh damping starts well above the boundary layer.
	1085	rdh = MIN( MAX( zi_ribulk * 1.3_wp, 10.0_wp * dz(1) ), &
	1086	0.8_wp * zu(nzt), rayleigh_damping_height )
	1087	!
	1088	!-- Update Rayleigh damping factor
	1089	DO k = nzb+1, nzt
	1090	IF ( zu(k) >= rdh ) THEN
	1091	rdf(k) = rayleigh_damping_factor * &
	1092	( SIN( pi * 0.5_wp * ( zu(k) - rdh ) &
	1093	/ ( zu(nzt) - rdh ) ) &
	1094	)**2
	1095	ENDIF
	1096	ENDDO
	1097	rdf_sc = rdf
	1098
	1099	ENDIF
	1100
	1101	END SUBROUTINE adjust_sponge_layer
	1102
	1103	!------------------------------------------------------------------------------!
	1104	! Description:
	1105	! ------------
	1106	!> Performs consistency checks
	1107	!------------------------------------------------------------------------------!
	1108	SUBROUTINE nesting_offl_check_parameters
	1109
	1110	IMPLICIT NONE
	1111	!
[4169]	1112	!-- Check if offline nesting is applied in nested child domain.
[3579]	1113	IF ( nesting_offline .AND. child_domain ) THEN
	1114	message_string = 'Offline nesting is only applicable in root model.'
[4169]	1115	CALL message( 'offline_nesting_check_parameters', 'PA0622', 1, 2, 0, 6, 0 )
	1116	ENDIF
[3347]	1117
	1118	END SUBROUTINE nesting_offl_check_parameters
	1119
	1120	!------------------------------------------------------------------------------!
	1121	! Description:
	1122	! ------------
	1123	!> Reads the parameter list nesting_offl_parameters
	1124	!------------------------------------------------------------------------------!
	1125	SUBROUTINE nesting_offl_parin
	1126
	1127	IMPLICIT NONE
	1128
	1129	CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file
	1130
	1131
	1132	NAMELIST /nesting_offl_parameters/ nesting_offline
	1133
	1134	line = ' '
	1135
	1136	!
	1137	!-- Try to find stg package
	1138	REWIND ( 11 )
	1139	line = ' '
	1140	DO WHILE ( INDEX( line, '&nesting_offl_parameters' ) == 0 )
	1141	READ ( 11, '(A)', END=20 ) line
	1142	ENDDO
	1143	BACKSPACE ( 11 )
	1144
	1145	!
	1146	!-- Read namelist
	1147	READ ( 11, nesting_offl_parameters, ERR = 10, END = 20 )
	1148
	1149	GOTO 20
	1150
	1151	10 BACKSPACE( 11 )
	1152	READ( 11 , '(A)') line
	1153	CALL parin_fail_message( 'nesting_offl_parameters', line )
	1154
	1155	20 CONTINUE
	1156
	1157
	1158	END SUBROUTINE nesting_offl_parin
	1159
	1160	!------------------------------------------------------------------------------!
	1161	! Description:
	1162	! ------------
	1163	!> Writes information about offline nesting into HEADER file
	1164	!------------------------------------------------------------------------------!
	1165	SUBROUTINE nesting_offl_header ( io )
	1166
	1167	IMPLICIT NONE
	1168
	1169	INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file
	1170
	1171	WRITE ( io, 1 )
	1172	IF ( nesting_offline ) THEN
	1173	WRITE ( io, 3 )
	1174	ELSE
	1175	WRITE ( io, 2 )
	1176	ENDIF
	1177
	1178	1 FORMAT (//' Offline nesting in COSMO model:'/ &
	1179	' -------------------------------'/)
	1180	2 FORMAT (' --> No offlince nesting is used (default) ')
	1181	3 FORMAT (' --> Offlince nesting is used. Boundary data is read from dynamic input file ')
	1182
	1183	END SUBROUTINE nesting_offl_header
	1184
	1185	!------------------------------------------------------------------------------!
	1186	! Description:
	1187	! ------------
	1188	!> Allocate arrays used to read boundary data from NetCDF file and initialize
	1189	!> boundary data.
	1190	!------------------------------------------------------------------------------!
	1191	SUBROUTINE nesting_offl_init
	1192
	1193	USE netcdf_data_input_mod, &
	1194	ONLY: netcdf_data_input_offline_nesting
	1195
	1196	IMPLICIT NONE
[3737]	1197
	1198	INTEGER(iwp) :: n !< running index for chemical species
[3347]	1199
	1200
	1201	!-- Allocate arrays for geostrophic wind components. Arrays will
[3404]	1202	!-- incorporate 2 time levels in order to interpolate in between.
	1203	ALLOCATE( nest_offl%ug(0:1,1:nzt) )
	1204	ALLOCATE( nest_offl%vg(0:1,1:nzt) )
[3347]	1205	!
[4125]	1206	!-- Allocate arrays for reading left/right boundary values. Arrays will
	1207	!-- incorporate 2 time levels in order to interpolate in between. If the core has
	1208	!-- no boundary, allocate a dummy array, in order to enable netcdf parallel
	1209	!-- access. Dummy arrays will be allocated with dimension length zero.
[3347]	1210	IF ( bc_dirichlet_l ) THEN
	1211	ALLOCATE( nest_offl%u_left(0:1,nzb+1:nzt,nys:nyn) )
	1212	ALLOCATE( nest_offl%v_left(0:1,nzb+1:nzt,nysv:nyn) )
	1213	ALLOCATE( nest_offl%w_left(0:1,nzb+1:nzt-1,nys:nyn) )
	1214	IF ( humidity ) ALLOCATE( nest_offl%q_left(0:1,nzb+1:nzt,nys:nyn) )
	1215	IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_left(0:1,nzb+1:nzt,nys:nyn) )
[3737]	1216	IF ( air_chemistry ) ALLOCATE( nest_offl%chem_left(0:1,nzb+1:nzt,nys:nyn,&
	1217	1:UBOUND( chem_species, 1 )) )
[4125]	1218	ELSE
	1219	ALLOCATE( nest_offl%u_left(1:1,1:1,1:1) )
	1220	ALLOCATE( nest_offl%v_left(1:1,1:1,1:1) )
	1221	ALLOCATE( nest_offl%w_left(1:1,1:1,1:1) )
	1222	IF ( humidity ) ALLOCATE( nest_offl%q_left(1:1,1:1,1:1) )
	1223	IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_left(1:1,1:1,1:1) )
	1224	IF ( air_chemistry ) ALLOCATE( nest_offl%chem_left(1:1,1:1,1:1, &
	1225	1:UBOUND( chem_species, 1 )) )
[3347]	1226	ENDIF
	1227	IF ( bc_dirichlet_r ) THEN
	1228	ALLOCATE( nest_offl%u_right(0:1,nzb+1:nzt,nys:nyn) )
	1229	ALLOCATE( nest_offl%v_right(0:1,nzb+1:nzt,nysv:nyn) )
	1230	ALLOCATE( nest_offl%w_right(0:1,nzb+1:nzt-1,nys:nyn) )
	1231	IF ( humidity ) ALLOCATE( nest_offl%q_right(0:1,nzb+1:nzt,nys:nyn) )
	1232	IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_right(0:1,nzb+1:nzt,nys:nyn) )
[3737]	1233	IF ( air_chemistry ) ALLOCATE( nest_offl%chem_right(0:1,nzb+1:nzt,nys:nyn,&
	1234	1:UBOUND( chem_species, 1 )) )
[4125]	1235	ELSE
	1236	ALLOCATE( nest_offl%u_right(1:1,1:1,1:1) )
	1237	ALLOCATE( nest_offl%v_right(1:1,1:1,1:1) )
	1238	ALLOCATE( nest_offl%w_right(1:1,1:1,1:1) )
	1239	IF ( humidity ) ALLOCATE( nest_offl%q_right(1:1,1:1,1:1) )
	1240	IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_right(1:1,1:1,1:1) )
	1241	IF ( air_chemistry ) ALLOCATE( nest_offl%chem_right(1:1,1:1,1:1, &
	1242	1:UBOUND( chem_species, 1 )) )
[3347]	1243	ENDIF
[4125]	1244	!
	1245	!-- Allocate arrays for reading north/south boundary values. Arrays will
	1246	!-- incorporate 2 time levels in order to interpolate in between. If the core has
	1247	!-- no boundary, allocate a dummy array, in order to enable netcdf parallel
	1248	!-- access. Dummy arrays will be allocated with dimension length zero.
[3347]	1249	IF ( bc_dirichlet_n ) THEN
	1250	ALLOCATE( nest_offl%u_north(0:1,nzb+1:nzt,nxlu:nxr) )
	1251	ALLOCATE( nest_offl%v_north(0:1,nzb+1:nzt,nxl:nxr) )
	1252	ALLOCATE( nest_offl%w_north(0:1,nzb+1:nzt-1,nxl:nxr) )
	1253	IF ( humidity ) ALLOCATE( nest_offl%q_north(0:1,nzb+1:nzt,nxl:nxr) )
	1254	IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_north(0:1,nzb+1:nzt,nxl:nxr) )
[3737]	1255	IF ( air_chemistry ) ALLOCATE( nest_offl%chem_north(0:1,nzb+1:nzt,nxl:nxr,&
	1256	1:UBOUND( chem_species, 1 )) )
[4125]	1257	ELSE
	1258	ALLOCATE( nest_offl%u_north(1:1,1:1,1:1) )
	1259	ALLOCATE( nest_offl%v_north(1:1,1:1,1:1) )
	1260	ALLOCATE( nest_offl%w_north(1:1,1:1,1:1) )
	1261	IF ( humidity ) ALLOCATE( nest_offl%q_north(1:1,1:1,1:1) )
	1262	IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_north(1:1,1:1,1:1) )
	1263	IF ( air_chemistry ) ALLOCATE( nest_offl%chem_north(1:1,1:1,1:1, &
	1264	1:UBOUND( chem_species, 1 )) )
[3347]	1265	ENDIF
	1266	IF ( bc_dirichlet_s ) THEN
	1267	ALLOCATE( nest_offl%u_south(0:1,nzb+1:nzt,nxlu:nxr) )
	1268	ALLOCATE( nest_offl%v_south(0:1,nzb+1:nzt,nxl:nxr) )
	1269	ALLOCATE( nest_offl%w_south(0:1,nzb+1:nzt-1,nxl:nxr) )
	1270	IF ( humidity ) ALLOCATE( nest_offl%q_south(0:1,nzb+1:nzt,nxl:nxr) )
	1271	IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_south(0:1,nzb+1:nzt,nxl:nxr) )
[3737]	1272	IF ( air_chemistry ) ALLOCATE( nest_offl%chem_south(0:1,nzb+1:nzt,nxl:nxr,&
	1273	1:UBOUND( chem_species, 1 )) )
[4125]	1274	ELSE
	1275	ALLOCATE( nest_offl%u_south(1:1,1:1,1:1) )
	1276	ALLOCATE( nest_offl%v_south(1:1,1:1,1:1) )
	1277	ALLOCATE( nest_offl%w_south(1:1,1:1,1:1) )
	1278	IF ( humidity ) ALLOCATE( nest_offl%q_south(1:1,1:1,1:1) )
	1279	IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_south(1:1,1:1,1:1) )
	1280	IF ( air_chemistry ) ALLOCATE( nest_offl%chem_south(1:1,1:1,1:1, &
	1281	1:UBOUND( chem_species, 1 )) )
[3347]	1282	ENDIF
[4125]	1283	!
	1284	!-- Allocate arrays for reading data at the top boundary. In contrast to the
	1285	!-- lateral boundaries, every core reads these data so that no dummy
	1286	!-- arrays need to be allocated.
[3347]	1287	ALLOCATE( nest_offl%u_top(0:1,nys:nyn,nxlu:nxr) )
	1288	ALLOCATE( nest_offl%v_top(0:1,nysv:nyn,nxl:nxr) )
	1289	ALLOCATE( nest_offl%w_top(0:1,nys:nyn,nxl:nxr) )
	1290	IF ( humidity ) ALLOCATE( nest_offl%q_top(0:1,nys:nyn,nxl:nxr) )
	1291	IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_top(0:1,nys:nyn,nxl:nxr) )
[3737]	1292	IF ( air_chemistry ) ALLOCATE( nest_offl%chem_top(0:1,nys:nyn,nxl:nxr, &
	1293	1:UBOUND( chem_species, 1 )) )
[3347]	1294	!
[3737]	1295	!-- For chemical species, create the names of the variables. This is necessary
	1296	!-- to identify the respective variable and write it onto the correct array
	1297	!-- in the chem_species datatype.
	1298	IF ( air_chemistry ) THEN
	1299	ALLOCATE( nest_offl%chem_from_file_l(1:UBOUND( chem_species, 1 )) )
	1300	ALLOCATE( nest_offl%chem_from_file_n(1:UBOUND( chem_species, 1 )) )
	1301	ALLOCATE( nest_offl%chem_from_file_r(1:UBOUND( chem_species, 1 )) )
	1302	ALLOCATE( nest_offl%chem_from_file_s(1:UBOUND( chem_species, 1 )) )
	1303	ALLOCATE( nest_offl%chem_from_file_t(1:UBOUND( chem_species, 1 )) )
	1304
	1305	ALLOCATE( nest_offl%var_names_chem_l(1:UBOUND( chem_species, 1 )) )
	1306	ALLOCATE( nest_offl%var_names_chem_n(1:UBOUND( chem_species, 1 )) )
	1307	ALLOCATE( nest_offl%var_names_chem_r(1:UBOUND( chem_species, 1 )) )
	1308	ALLOCATE( nest_offl%var_names_chem_s(1:UBOUND( chem_species, 1 )) )
	1309	ALLOCATE( nest_offl%var_names_chem_t(1:UBOUND( chem_species, 1 )) )
	1310	!
	1311	!-- Initialize flags that indicate whether the variable is on file or
	1312	!-- not. Please note, this is only necessary for chemistry variables.
	1313	nest_offl%chem_from_file_l(:) = .FALSE.
	1314	nest_offl%chem_from_file_n(:) = .FALSE.
	1315	nest_offl%chem_from_file_r(:) = .FALSE.
	1316	nest_offl%chem_from_file_s(:) = .FALSE.
	1317	nest_offl%chem_from_file_t(:) = .FALSE.
	1318
	1319	DO n = 1, UBOUND( chem_species, 1 )
	1320	nest_offl%var_names_chem_l(n) = nest_offl%char_l // &
	1321	TRIM(chem_species(n)%name)
	1322	nest_offl%var_names_chem_n(n) = nest_offl%char_n // &
	1323	TRIM(chem_species(n)%name)
	1324	nest_offl%var_names_chem_r(n) = nest_offl%char_r // &
	1325	TRIM(chem_species(n)%name)
	1326	nest_offl%var_names_chem_s(n) = nest_offl%char_s // &
	1327	TRIM(chem_species(n)%name)
	1328	nest_offl%var_names_chem_t(n) = nest_offl%char_t // &
	1329	TRIM(chem_species(n)%name)
	1330	ENDDO
	1331	ENDIF
	1332	!
[3347]	1333	!-- Read COSMO data at lateral and top boundaries
	1334	CALL netcdf_data_input_offline_nesting
	1335	!
[4169]	1336	!-- Check if sufficient time steps are provided to cover the entire
	1337	!-- simulation. Note, dynamic input is only required for the 3D simulation,
	1338	!-- not for the soil/wall spinup. However, as the spinup time is added
	1339	!-- to the end_time, this must be considered here.
	1340	IF ( end_time - spinup_time > nest_offl%time(nest_offl%nt-1) ) THEN
	1341	message_string = 'end_time > provided time in offline nesting.'
	1342	CALL message( 'offline_nesting_check_parameters', 'PA0183', &
	1343	1, 2, 0, 6, 0 )
	1344	ENDIF
	1345	!
[3891]	1346	!-- Initialize boundary data. Please note, do not initialize boundaries in
	1347	!-- case of restart runs. This case the boundaries are already initialized
	1348	!-- and the boundary data from file would be on the wrong time level.
	1349	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
	1350	IF ( bc_dirichlet_l ) THEN
	1351	u(nzb+1:nzt,nys:nyn,0) = nest_offl%u_left(0,nzb+1:nzt,nys:nyn)
	1352	v(nzb+1:nzt,nysv:nyn,-1) = nest_offl%v_left(0,nzb+1:nzt,nysv:nyn)
	1353	w(nzb+1:nzt-1,nys:nyn,-1) = nest_offl%w_left(0,nzb+1:nzt-1,nys:nyn)
	1354	IF ( .NOT. neutral ) pt(nzb+1:nzt,nys:nyn,-1) = &
	1355	nest_offl%pt_left(0,nzb+1:nzt,nys:nyn)
	1356	IF ( humidity ) q(nzb+1:nzt,nys:nyn,-1) = &
	1357	nest_offl%q_left(0,nzb+1:nzt,nys:nyn)
	1358	IF ( air_chemistry ) THEN
	1359	DO n = 1, UBOUND( chem_species, 1 )
	1360	IF( nest_offl%chem_from_file_l(n) ) THEN
	1361	chem_species(n)%conc(nzb+1:nzt,nys:nyn,-1) = &
	1362	nest_offl%chem_left(0,nzb+1:nzt,nys:nyn,n)
	1363	ENDIF
	1364	ENDDO
	1365	ENDIF
[3737]	1366	ENDIF
[3891]	1367	IF ( bc_dirichlet_r ) THEN
	1368	u(nzb+1:nzt,nys:nyn,nxr+1) = nest_offl%u_right(0,nzb+1:nzt,nys:nyn)
	1369	v(nzb+1:nzt,nysv:nyn,nxr+1) = nest_offl%v_right(0,nzb+1:nzt,nysv:nyn)
	1370	w(nzb+1:nzt-1,nys:nyn,nxr+1) = nest_offl%w_right(0,nzb+1:nzt-1,nys:nyn)
	1371	IF ( .NOT. neutral ) pt(nzb+1:nzt,nys:nyn,nxr+1) = &
	1372	nest_offl%pt_right(0,nzb+1:nzt,nys:nyn)
	1373	IF ( humidity ) q(nzb+1:nzt,nys:nyn,nxr+1) = &
	1374	nest_offl%q_right(0,nzb+1:nzt,nys:nyn)
	1375	IF ( air_chemistry ) THEN
	1376	DO n = 1, UBOUND( chem_species, 1 )
	1377	IF( nest_offl%chem_from_file_r(n) ) THEN
	1378	chem_species(n)%conc(nzb+1:nzt,nys:nyn,nxr+1) = &
	1379	nest_offl%chem_right(0,nzb+1:nzt,nys:nyn,n)
	1380	ENDIF
	1381	ENDDO
	1382	ENDIF
[3737]	1383	ENDIF
[3891]	1384	IF ( bc_dirichlet_s ) THEN
	1385	u(nzb+1:nzt,-1,nxlu:nxr) = nest_offl%u_south(0,nzb+1:nzt,nxlu:nxr)
	1386	v(nzb+1:nzt,0,nxl:nxr) = nest_offl%v_south(0,nzb+1:nzt,nxl:nxr)
	1387	w(nzb+1:nzt-1,-1,nxl:nxr) = nest_offl%w_south(0,nzb+1:nzt-1,nxl:nxr)
	1388	IF ( .NOT. neutral ) pt(nzb+1:nzt,-1,nxl:nxr) = &
	1389	nest_offl%pt_south(0,nzb+1:nzt,nxl:nxr)
	1390	IF ( humidity ) q(nzb+1:nzt,-1,nxl:nxr) = &
	1391	nest_offl%q_south(0,nzb+1:nzt,nxl:nxr)
	1392	IF ( air_chemistry ) THEN
	1393	DO n = 1, UBOUND( chem_species, 1 )
	1394	IF( nest_offl%chem_from_file_s(n) ) THEN
	1395	chem_species(n)%conc(nzb+1:nzt,-1,nxl:nxr) = &
	1396	nest_offl%chem_south(0,nzb+1:nzt,nxl:nxr,n)
	1397	ENDIF
	1398	ENDDO
	1399	ENDIF
[3737]	1400	ENDIF
[3891]	1401	IF ( bc_dirichlet_n ) THEN
	1402	u(nzb+1:nzt,nyn+1,nxlu:nxr) = nest_offl%u_north(0,nzb+1:nzt,nxlu:nxr)
	1403	v(nzb+1:nzt,nyn+1,nxl:nxr) = nest_offl%v_north(0,nzb+1:nzt,nxl:nxr)
	1404	w(nzb+1:nzt-1,nyn+1,nxl:nxr) = nest_offl%w_north(0,nzb+1:nzt-1,nxl:nxr)
	1405	IF ( .NOT. neutral ) pt(nzb+1:nzt,nyn+1,nxl:nxr) = &
	1406	nest_offl%pt_north(0,nzb+1:nzt,nxl:nxr)
	1407	IF ( humidity ) q(nzb+1:nzt,nyn+1,nxl:nxr) = &
	1408	nest_offl%q_north(0,nzb+1:nzt,nxl:nxr)
	1409	IF ( air_chemistry ) THEN
	1410	DO n = 1, UBOUND( chem_species, 1 )
	1411	IF( nest_offl%chem_from_file_n(n) ) THEN
	1412	chem_species(n)%conc(nzb+1:nzt,nyn+1,nxl:nxr) = &
	1413	nest_offl%chem_north(0,nzb+1:nzt,nxl:nxr,n)
	1414	ENDIF
	1415	ENDDO
	1416	ENDIF
[3737]	1417	ENDIF
[3891]	1418	!
	1419	!-- Initialize geostrophic wind components. Actually this is already done in
	1420	!-- init_3d_model when initializing_action = 'inifor', however, in speical
	1421	!-- case of user-defined initialization this will be done here again, in
	1422	!-- order to have a consistent initialization.
	1423	ug(nzb+1:nzt) = nest_offl%ug(0,nzb+1:nzt)
	1424	vg(nzb+1:nzt) = nest_offl%vg(0,nzb+1:nzt)
	1425	!
	1426	!-- Set bottom and top boundary condition for geostrophic wind components
	1427	ug(nzt+1) = ug(nzt)
	1428	vg(nzt+1) = vg(nzt)
	1429	ug(nzb) = ug(nzb+1)
	1430	vg(nzb) = vg(nzb+1)
	1431	ENDIF
[3347]	1432	!
	1433	!-- After boundary data is initialized, mask topography at the
	1434	!-- boundaries for the velocity components.
	1435	u = MERGE( u, 0.0_wp, BTEST( wall_flags_0, 1 ) )
	1436	v = MERGE( v, 0.0_wp, BTEST( wall_flags_0, 2 ) )
	1437	w = MERGE( w, 0.0_wp, BTEST( wall_flags_0, 3 ) )
[3891]	1438	!
	1439	!-- Initial calculation of the boundary layer depth from the prescribed
	1440	!-- boundary data. This is requiered for initialize the synthetic turbulence
	1441	!-- generator correctly.
[4022]	1442	CALL nesting_offl_calc_zi
[3347]	1443
	1444	!
[3891]	1445	!-- After boundary data is initialized, ensure mass conservation. Not
	1446	!-- necessary in restart runs.
	1447	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
	1448	CALL nesting_offl_mass_conservation
	1449	ENDIF
[3347]	1450
	1451	END SUBROUTINE nesting_offl_init
	1452
	1453	!------------------------------------------------------------------------------!
	1454	! Description:
	1455	!------------------------------------------------------------------------------!
	1456	!> Interpolation function, used to interpolate boundary data in time.
	1457	!------------------------------------------------------------------------------!
	1458	FUNCTION interpolate_in_time( var_t1, var_t2, fac )
	1459
	1460	USE kinds
	1461
	1462	IMPLICIT NONE
	1463
	1464	REAL(wp) :: interpolate_in_time !< time-interpolated boundary value
	1465	REAL(wp) :: var_t1 !< boundary value at t1
	1466	REAL(wp) :: var_t2 !< boundary value at t2
	1467	REAL(wp) :: fac !< interpolation factor
	1468
	1469	interpolate_in_time = ( 1.0_wp - fac ) * var_t1 + fac * var_t2
	1470
	1471	END FUNCTION interpolate_in_time
	1472
	1473
	1474
	1475	END MODULE nesting_offl_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |