Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |