Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |