Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |