Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |