Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |