Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |