Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |