Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |