Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |