Home

Context Navigation

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

Last change on this file since 4102 was 4079, checked in by suehring, 5 years ago
Implementation of a monotonic flux limiter for vertical advection term in Wicker-Skamarock scheme. The flux limiter is currently only applied for passive scalars (passive scalar, chemical species, aerosols) within the region up to the highest topography, in order to avoid the built-up of large concentrations within poorly resolved cavities in urban environments. To enable the limiter monotonic_limiter_z = .T. must be set. Note, the limiter is currently only implemented for the cache-optimized version of advec_ws. Further changes in offline nesting: Set boundary condition for w at nzt+1 at all lateral boundaries (even though these won't enter the numerical solution), in order to avoid high vertical velocities in the run-control file which might built-up due to the mass-conservation; bugfix in offline nesting for chemical species
Property svn:keywords set to `Id`
File size: 61.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |