Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |