Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |