Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |