Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |