Home

Context Navigation

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

Last change on this file since 3890 was 3885, checked in by kanani, 5 years ago
restructure/add location/debug messages
Property svn:keywords set to `Id`
File size: 56.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |