Home

Context Navigation

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

Last change on this file since 4178 was 4169, checked in by suehring, 5 years ago
Check added
Property svn:keywords set to `Id`
File size: 66.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |