Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |