Home

Context Navigation

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

Last change on this file since 3704 was 3704, checked in by suehring, 5 years ago
Revision of virtual-measurement module and data output enabled. Further, post-processing tool added to merge distributed virtually sampled data and to output it into NetCDF files.
Property svn:keywords set to `Id`
File size: 47.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |