Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |