Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |