Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |