Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |