Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |