Home

Context Navigation

source: palm/trunk/SOURCE/large_scale_forcing_nudging_mod.f90 @ 2956

Last change on this file since 2956 was 2938, checked in by suehring, 7 years ago
Nesting in RANS-LES and RANS-RANS mode enabled; synthetic turbulence generator at all lateral boundaries in nesting or non-cyclic forcing mode; revised Inifor initialization in nesting mode
Property svn:keywords set to `Id`
File size: 89.1 KB

Line
1	!> @file large_scale_forcing_nudging_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-2018 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! ------------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: large_scale_forcing_nudging_mod.f90 2938 2018-03-27 15:52:42Z Giersch $
27	! Further improvements for nesting in larger-scale model
28	!
29	! 2863 2018-03-08 11:36:25Z suehring
30	! Corrected "Former revisions" section
31	!
32	! 2696 2017-12-14 17:12:51Z kanani
33	! Change in file header (GPL part)
34	! Forcing with larger-scale models implemented (MS)
35	!
36	! 2342 2017-08-08 11:00:43Z boeske
37	! fixed check if surface forcing data is available until end of simulation
38	!
39	! 2320 2017-07-21 12:47:43Z suehring
40	! initial revision
41	!
42	! Description:
43	! ------------
44	!> Calculates large scale forcings (geostrophic wind and subsidence velocity) as
45	!> well as surfaces fluxes dependent on time given in an external file (LSF_DATA).
46	!> Moreover, module contains nudging routines, where u, v, pt and q are nudged
47	!> to given profiles on a relaxation timescale tnudge.
48	!> Profiles are read in from NUDGING_DATA.
49	!> Code is based on Neggers et al. (2012) and also in parts on DALES and UCLA-LES.
50	!> @todo: Revise reading of ASCII-files
51	!> @todo: Remove unused variables and control flags
52	!> @todo: Revise large-scale facing of surface variables
53	!> @todo: Revise control flags lsf_exception, lsf_surf, lsf_vert, etc.
54	!--------------------------------------------------------------------------------!
55	MODULE lsf_nudging_mod
56
57	USE arrays_3d, &
58	ONLY: dzw, e, heatflux_input_conversion, pt, pt_init, q, q_init, s, &
59	tend, u, u_init, ug, v, v_init, vg, w, w_subs, &
60	waterflux_input_conversion, zu, zw
61
62	USE control_parameters, &
63	ONLY: bc_lr, bc_ns, bc_pt_b, bc_q_b, constant_diffusion, &
64	data_output_pr, dt_3d, end_time, forcing, &
65	force_bound_l, force_bound_n, force_bound_r, force_bound_s, &
66	humidity, intermediate_timestep_count, ibc_pt_b, ibc_q_b, &
67	large_scale_forcing, large_scale_subsidence, lsf_surf, lsf_vert,&
68	lsf_exception, message_string, neutral, nudging, passive_scalar,&
69	pt_surface, ocean, q_surface, surface_pressure, topography, &
70	use_subsidence_tendencies
71
72	USE grid_variables
73
74	USE pegrid
75
76	USE indices, &
77	ONLY: nbgp, ngp_sums_ls, nx, nxl, nxlg, nxlu, nxr, nxrg, ny, nys, &
78	nysv, nysg, nyn, nyng, nzb, nz, nzt, wall_flags_0
79
80	USE kinds
81
82	USE surface_mod, &
83	ONLY: surf_def_h, surf_lsm_h, surf_usm_h
84
85	USE statistics, &
86	ONLY: hom, statistic_regions, sums_ls_l, weight_substep
87
88	USE netcdf_data_input_mod, &
89	ONLY: force, netcdf_data_input_interpolate
90
91	INTEGER(iwp) :: nlsf = 1000 !< maximum number of profiles in LSF_DATA (large scale forcing)
92	INTEGER(iwp) :: ntnudge = 1000 !< maximum number of profiles in NUDGING_DATA (nudging)
93
94	REAL(wp) :: d_area_t
95
96	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ptnudge !< vertical profile of pot. temperature interpolated to vertical grid (nudging)
97	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: qnudge !< vertical profile of specific humidity interpolated to vertical grid (nudging)
98	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tnudge !< vertical profile of nudging time scale interpolated to vertical grid (nudging)
99	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: td_lsa_lpt !< temperature tendency due to large scale advection (large scale forcing)
100	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: td_lsa_q !< specific humidity tendency due to large scale advection (large scale forcing)
101	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: td_sub_lpt !< temperature tendency due to subsidence/ascent (large scale forcing)
102	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: td_sub_q !< specific humidity tendency due to subsidence/ascent (large scale forcing)
103	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ug_vert !< vertical profile of geostrophic wind component in x-direction interpolated to vertical grid (large scale forcing)
104	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: unudge !< vertical profile of wind component in x-direction interpolated to vertical grid (nudging)
105	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: vnudge !< vertical profile of wind component in y-direction interpolated to vertical grid (nudging)
106	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: vg_vert !< vertical profile of geostrophic wind component in y-direction interpolated to vertical grid (large scale forcing)
107	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: wnudge !< vertical profile of subsidence/ascent velocity interpolated to vertical grid (nudging) ???
108	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: wsubs_vert !< vertical profile of wind component in z-direction interpolated to vertical grid (nudging) ???
109
110	REAL(wp), DIMENSION(:), ALLOCATABLE :: shf_surf !< time-dependent surface sensible heat flux (large scale forcing)
111	REAL(wp), DIMENSION(:), ALLOCATABLE :: timenudge !< times at which vertical profiles are defined in NUDGING_DATA (nudging)
112	REAL(wp), DIMENSION(:), ALLOCATABLE :: time_surf !< times at which surface values/fluxes are defined in LSF_DATA (large scale forcing)
113	REAL(wp), DIMENSION(:), ALLOCATABLE :: time_vert !< times at which vertical profiles are defined in LSF_DATA (large scale forcing)
114
115	REAL(wp), DIMENSION(:), ALLOCATABLE :: tmp_tnudge !< current nudging time scale
116
117	REAL(wp), DIMENSION(:), ALLOCATABLE :: p_surf !< time-dependent surface pressure (large scale forcing)
118	REAL(wp), DIMENSION(:), ALLOCATABLE :: pt_surf !< time-dependent surface temperature (large scale forcing)
119	REAL(wp), DIMENSION(:), ALLOCATABLE :: qsws_surf !< time-dependent surface latent heat flux (large scale forcing)
120	REAL(wp), DIMENSION(:), ALLOCATABLE :: q_surf !< time-dependent surface specific humidity (large scale forcing)
121
122	SAVE
123	PRIVATE
124	!
125	!-- Public subroutines
126	PUBLIC ls_forcing_surf, ls_forcing_vert, ls_advec, lsf_init, &
127	lsf_nudging_check_parameters, nudge_init, &
128	lsf_nudging_check_data_output_pr, lsf_nudging_header, &
129	calc_tnudge, nudge, nudge_ref, forcing_bc_mass_conservation, &
130	forcing_bc
131	!
132	!-- Public variables
133	PUBLIC qsws_surf, shf_surf, td_lsa_lpt, td_lsa_q, td_sub_lpt, &
134	td_sub_q, time_vert, force
135
136
137	INTERFACE ls_advec
138	MODULE PROCEDURE ls_advec
139	MODULE PROCEDURE ls_advec_ij
140	END INTERFACE ls_advec
141
142	INTERFACE nudge
143	MODULE PROCEDURE nudge
144	MODULE PROCEDURE nudge_ij
145	END INTERFACE nudge
146
147	CONTAINS
148
149
150	!------------------------------------------------------------------------------!
151	! Description:
152	! ------------
153	!> @todo Missing subroutine description.
154	!------------------------------------------------------------------------------!
155	SUBROUTINE forcing_bc_mass_conservation
156
157	USE control_parameters, &
158	ONLY: volume_flow
159
160	IMPLICIT NONE
161
162	INTEGER(iwp) :: i !<
163	INTEGER(iwp) :: j !<
164	INTEGER(iwp) :: k !<
165
166	REAL(wp) :: w_correct !<
167	REAL(wp), DIMENSION(1:3) :: volume_flow_l !<
168
169	volume_flow = 0.0_wp
170	volume_flow_l = 0.0_wp
171
172	d_area_t = 1.0_wp / ( ( nx + 1 ) * dx * ( ny + 1 ) * dy )
173
174	IF ( force_bound_l ) THEN
175	i = nxl
176	DO j = nys, nyn
177	DO k = nzb+1, nzt
178	volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k) * dy &
179	* MERGE( 1.0_wp, 0.0_wp, &
180	BTEST( wall_flags_0(k,j,i), 1 ) )
181	ENDDO
182	ENDDO
183	ENDIF
184	IF ( force_bound_r ) THEN
185	i = nxr+1
186	DO j = nys, nyn
187	DO k = nzb+1, nzt
188	volume_flow_l(1) = volume_flow_l(1) - u(k,j,i) * dzw(k) * dy &
189	* MERGE( 1.0_wp, 0.0_wp, &
190	BTEST( wall_flags_0(k,j,i), 1 ) )
191	ENDDO
192	ENDDO
193	ENDIF
194	IF ( force_bound_s ) THEN
195	j = nys
196	DO i = nxl, nxr
197	DO k = nzb+1, nzt
198	volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k) * dx &
199	* MERGE( 1.0_wp, 0.0_wp, &
200	BTEST( wall_flags_0(k,j,i), 2 ) )
201	ENDDO
202	ENDDO
203	ENDIF
204	IF ( force_bound_n ) THEN
205	j = nyn+1
206	DO i = nxl, nxr
207	DO k = nzb+1, nzt
208	volume_flow_l(2) = volume_flow_l(2) - v(k,j,i) * dzw(k) * dx &
209	* MERGE( 1.0_wp, 0.0_wp, &
210	BTEST( wall_flags_0(k,j,i), 2 ) )
211	ENDDO
212	ENDDO
213	ENDIF
214	!
215	!-- Top boundary
216	k = nzt
217	DO i = nxl, nxr
218	DO j = nys, nyn
219	volume_flow_l(3) = volume_flow_l(3) - w(k,j,i) * dx * dy
220	ENDDO
221	ENDDO
222
223	#if defined( __parallel )
224	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
225	CALL MPI_ALLREDUCE( volume_flow_l, volume_flow, 3, MPI_REAL, MPI_SUM, &
226	comm2d, ierr )
227	#else
228	volume_flow = volume_flow_l
229	#endif
230
231	w_correct = SUM( volume_flow ) * d_area_t
232
233	DO i = nxl, nxr
234	DO j = nys, nyn
235	DO k = nzt, nzt + 1
236	w(k,j,i) = w(k,j,i) + w_correct
237	ENDDO
238	ENDDO
239	ENDDO
240
241	write(9,*) "w correction", w_correct
242	flush(9)
243
244	END SUBROUTINE forcing_bc_mass_conservation
245
246
247	!------------------------------------------------------------------------------!
248	! Description:
249	! ------------
250	!> @todo Missing subroutine description.
251	!------------------------------------------------------------------------------!
252	SUBROUTINE forcing_bc
253
254	USE control_parameters, &
255	ONLY: force_bound_l, force_bound_n, force_bound_r, force_bound_s, &
256	humidity, neutral, passive_scalar, simulated_time
257
258	USE netcdf_data_input_mod, &
259	ONLY: force
260
261	IMPLICIT NONE
262
263	INTEGER(iwp) :: i !< running index x-direction
264	INTEGER(iwp) :: j !< running index y-direction
265	INTEGER(iwp) :: k !< running index z-direction
266	INTEGER(iwp) :: t !< running index for time levels
267
268	REAL(wp) :: ddt_lsf !< inverse value of time resolution of forcing data
269	REAL(wp) :: t_ref !< time past since last reference step
270
271	!
272	!-- If required, interpolate and/or extrapolate data vertically. This is
273	!-- required as Inifor outputs only equidistant vertical data.
274	IF ( ANY( zu(1:nzt+1) /= force%zu_atmos(1:force%nzu) ) ) THEN
275	IF ( .NOT. force%interpolated ) THEN
276
277	DO t = 0, 1
278	IF ( force_bound_l ) THEN
279	CALL netcdf_data_input_interpolate( force%u_left(t,:,:), &
280	zu(nzb+1:nzt+1), &
281	force%zu_atmos )
282	CALL netcdf_data_input_interpolate( force%v_left(t,:,:), &
283	zu(nzb+1:nzt+1), &
284	force%zu_atmos )
285	CALL netcdf_data_input_interpolate( force%w_left(t,:,:), &
286	zw(nzb+1:nzt+1), &
287	force%zw_atmos )
288	IF ( .NOT. neutral ) &
289	CALL netcdf_data_input_interpolate( force%pt_left(t,:,:),&
290	zu(nzb+1:nzt+1), &
291	force%zu_atmos )
292	IF ( humidity ) &
293	CALL netcdf_data_input_interpolate( force%q_left(t,:,:), &
294	zu(nzb+1:nzt+1), &
295	force%zu_atmos )
296	ENDIF
297	IF ( force_bound_r ) THEN
298	CALL netcdf_data_input_interpolate( force%u_right(t,:,:), &
299	zu(nzb+1:nzt+1), &
300	force%zu_atmos )
301	CALL netcdf_data_input_interpolate( force%v_right(t,:,:), &
302	zu(nzb+1:nzt+1), &
303	force%zu_atmos )
304	CALL netcdf_data_input_interpolate( force%w_right(t,:,:), &
305	zw(nzb+1:nzt+1), &
306	force%zw_atmos )
307	IF ( .NOT. neutral ) &
308	CALL netcdf_data_input_interpolate( force%pt_right(t,:,:),&
309	zu(nzb+1:nzt+1), &
310	force%zu_atmos )
311	IF ( humidity ) &
312	CALL netcdf_data_input_interpolate( force%q_right(t,:,:),&
313	zu(nzb+1:nzt+1), &
314	force%zu_atmos )
315	ENDIF
316	IF ( force_bound_n ) THEN
317	CALL netcdf_data_input_interpolate( force%u_north(t,:,:), &
318	zu(nzb+1:nzt+1), &
319	force%zu_atmos )
320	CALL netcdf_data_input_interpolate( force%v_north(t,:,:), &
321	zu(nzb+1:nzt+1), &
322	force%zu_atmos )
323	CALL netcdf_data_input_interpolate( force%w_north(t,:,:), &
324	zw(nzb+1:nzt+1), &
325	force%zw_atmos )
326	IF ( .NOT. neutral ) &
327	CALL netcdf_data_input_interpolate( force%pt_north(t,:,:),&
328	zu(nzb+1:nzt+1), &
329	force%zu_atmos )
330	IF ( humidity ) &
331	CALL netcdf_data_input_interpolate( force%q_north(t,:,:),&
332	zu(nzb+1:nzt+1), &
333	force%zu_atmos )
334	ENDIF
335	IF ( force_bound_s ) THEN
336	CALL netcdf_data_input_interpolate( force%u_south(t,:,:), &
337	zu(nzb+1:nzt+1), &
338	force%zu_atmos )
339	CALL netcdf_data_input_interpolate( force%v_south(t,:,:), &
340	zu(nzb+1:nzt+1), &
341	force%zu_atmos )
342	CALL netcdf_data_input_interpolate( force%w_south(t,:,:), &
343	zw(nzb+1:nzt+1), &
344	force%zw_atmos )
345	IF ( .NOT. neutral ) &
346	CALL netcdf_data_input_interpolate( force%pt_south(t,:,:),&
347	zu(nzb+1:nzt+1), &
348	force%zu_atmos )
349	IF ( humidity ) &
350	CALL netcdf_data_input_interpolate( force%q_south(t,:,:),&
351	zu(nzb+1:nzt+1), &
352	force%zu_atmos )
353	ENDIF
354	ENDDO
355	!
356	!-- Note, no interpolation of top boundary. Just use initial value.
357	!-- No physical meaningful extrapolation possible if only one layer is
358	!-- given.
359
360	force%interpolated = .TRUE.
361	ENDIF
362	ENDIF
363
364	!
365	!-- Calculate time interval of forcing data
366	ddt_lsf = 1.0_wp / ( force%time(force%tind_p) - force%time(force%tind) )
367	!
368	!-- Calculate reziproke time past since last reference step. Please note,
369	!-- as simulated time is still not updated, the actual time here is
370	!-- simulated time + dt_3d
371	t_ref = simulated_time + dt_3d - force%time(force%tind)
372
373	IF ( force_bound_l ) THEN
374
375	DO j = nys, nyn
376	DO k = nzb+1, nzt+1
377	u(k,j,nxlg:nxl) = force%u_left(0,k,j) + ddt_lsf * t_ref * &
378	( force%u_left(1,k,j) - force%u_left(0,k,j) ) * &
379	MERGE( 1.0_wp, 0.0_wp, &
380	BTEST( wall_flags_0(k,j,nxlg:nxl), 1 ) )
381	ENDDO
382	ENDDO
383
384	DO j = nys, nyn
385	DO k = nzb+1, nzt
386	w(k,j,nxlg:nxl-1) = force%w_left(0,k,j) + ddt_lsf * t_ref * &
387	( force%w_left(1,k,j) - force%w_left(0,k,j) ) * &
388	MERGE( 1.0_wp, 0.0_wp, &
389	BTEST( wall_flags_0(k,j,nxlg:nxl-1), 3 ) )
390	ENDDO
391	ENDDO
392
393	DO j = nysv, nyn
394	DO k = nzb+1, nzt+1
395	v(k,j,nxlg:nxl-1) = force%v_left(0,k,j) + ddt_lsf * t_ref * &
396	( force%v_left(1,k,j) - force%v_left(0,k,j) ) * &
397	MERGE( 1.0_wp, 0.0_wp, &
398	BTEST( wall_flags_0(k,j,nxlg:nxl-1), 2 ) )
399	ENDDO
400	ENDDO
401
402	IF ( .NOT. neutral ) THEN
403	DO j = nys, nyn
404	DO k = nzb+1, nzt+1
405	pt(k,j,nxlg:nxl-1) = force%pt_left(0,k,j) + ddt_lsf * &
406	t_ref * &
407	( force%pt_left(1,k,j) - force%pt_left(0,k,j) )
408
409	ENDDO
410	ENDDO
411	ENDIF
412
413	IF ( humidity ) THEN
414	DO j = nys, nyn
415	DO k = nzb+1, nzt+1
416	q(k,j,nxlg:nxl-1) = force%q_left(0,k,j) + ddt_lsf * &
417	t_ref * &
418	( force%q_left(1,k,j) - force%q_left(0,k,j) )
419
420	ENDDO
421	ENDDO
422	ENDIF
423
424	ENDIF
425
426	IF ( force_bound_r ) THEN
427
428	DO j = nys, nyn
429	DO k = nzb+1, nzt+1
430	u(k,j,nxr+1:nxrg) = force%u_right(0,k,j) + ddt_lsf * t_ref * &
431	( force%u_right(1,k,j) - force%u_right(0,k,j) ) * &
432	MERGE( 1.0_wp, 0.0_wp, &
433	BTEST( wall_flags_0(k,j,nxr+1:nxrg), 1 ) )
434
435	ENDDO
436	ENDDO
437	DO j = nys, nyn
438	DO k = nzb+1, nzt
439	w(k,j,nxr+1:nxrg) = force%w_right(0,k,j) + ddt_lsf * t_ref * &
440	( force%w_right(1,k,j) - force%w_right(0,k,j) ) * &
441	MERGE( 1.0_wp, 0.0_wp, &
442	BTEST( wall_flags_0(k,j,nxr+1:nxrg), 3 ) )
443	ENDDO
444	ENDDO
445
446	DO j = nysv, nyn
447	DO k = nzb+1, nzt+1
448	v(k,j,nxr+1:nxrg) = force%v_right(0,k,j) + ddt_lsf * t_ref * &
449	( force%v_right(1,k,j) - force%v_right(0,k,j) ) * &
450	MERGE( 1.0_wp, 0.0_wp, &
451	BTEST( wall_flags_0(k,j,nxr+1:nxrg), 2 ) )
452	ENDDO
453	ENDDO
454
455	IF ( .NOT. neutral ) THEN
456	DO j = nys, nyn
457	DO k = nzb+1, nzt+1
458	pt(k,j,nxr+1:nxrg) = force%pt_right(0,k,j) + ddt_lsf * &
459	t_ref * &
460	( force%pt_right(1,k,j) - force%pt_right(0,k,j) )
461
462	ENDDO
463	ENDDO
464	ENDIF
465
466	IF ( humidity ) THEN
467	DO j = nys, nyn
468	DO k = nzb+1, nzt+1
469	q(k,j,nxr+1:nxrg) = force%q_right(0,k,j) + ddt_lsf * &
470	t_ref * &
471	( force%q_right(1,k,j) - force%q_right(0,k,j) )
472
473	ENDDO
474	ENDDO
475	ENDIF
476
477	ENDIF
478
479	IF ( force_bound_s ) THEN
480
481	DO i = nxl, nxr
482	DO k = nzb+1, nzt+1
483	v(k,nysg:nys,i) = force%v_south(0,k,i) + ddt_lsf * t_ref * &
484	( force%v_south(1,k,i) - force%v_south(0,k,i) ) * &
485	MERGE( 1.0_wp, 0.0_wp, &
486	BTEST( wall_flags_0(k,nysg:nys,i), 2 ) )
487	ENDDO
488	ENDDO
489
490	DO i = nxl, nxr
491	DO k = nzb+1, nzt
492	w(k,nysg:nys-1,i) = force%w_south(0,k,i) + ddt_lsf * t_ref * &
493	( force%w_south(1,k,i) - force%w_south(0,k,i) ) * &
494	MERGE( 1.0_wp, 0.0_wp, &
495	BTEST( wall_flags_0(k,nysg:nys-1,i), 3 ) )
496	ENDDO
497	ENDDO
498
499	DO i = nxlu, nxr
500	DO k = nzb+1, nzt+1
501	u(k,nysg:nys-1,i) = force%u_south(0,k,i) + ddt_lsf * t_ref * &
502	( force%u_south(1,k,i) - force%u_south(0,k,i) ) * &
503	MERGE( 1.0_wp, 0.0_wp, &
504	BTEST( wall_flags_0(k,nysg:nys-1,i), 1 ) )
505	ENDDO
506	ENDDO
507
508	IF ( .NOT. neutral ) THEN
509	DO i = nxl, nxr
510	DO k = nzb+1, nzt+1
511	pt(k,nysg:nys-1,i) = force%pt_south(0,k,i) + ddt_lsf * &
512	t_ref * &
513	( force%pt_south(1,k,i) - force%pt_south(0,k,i) )
514
515	ENDDO
516	ENDDO
517	ENDIF
518
519	IF ( humidity ) THEN
520	DO i = nxl, nxr
521	DO k = nzb+1, nzt+1
522	q(k,nysg:nys-1,i) = force%q_south(0,k,i) + ddt_lsf * &
523	t_ref * &
524	( force%q_south(1,k,i) - force%q_south(0,k,i) )
525
526	ENDDO
527	ENDDO
528	ENDIF
529
530	ENDIF
531
532	IF ( force_bound_n ) THEN
533
534	DO i = nxl, nxr
535	DO k = nzb+1, nzt+1
536	v(k,nyn+1:nyng,i) = force%v_north(0,k,i) + ddt_lsf * t_ref * &
537	( force%v_north(1,k,i) - force%v_north(0,k,i) ) * &
538	MERGE( 1.0_wp, 0.0_wp, &
539	BTEST( wall_flags_0(k,nyn+1:nyng,i), 2 ) )
540	ENDDO
541	ENDDO
542	DO i = nxl, nxr
543	DO k = nzb+1, nzt
544	w(k,nyn+1:nyng,i) = force%w_north(0,k,i) + ddt_lsf * t_ref * &
545	( force%w_north(1,k,i) - force%w_north(0,k,i) ) * &
546	MERGE( 1.0_wp, 0.0_wp, &
547	BTEST( wall_flags_0(k,nyn+1:nyng,i), 3 ) )
548	ENDDO
549	ENDDO
550
551	DO i = nxlu, nxr
552	DO k = nzb+1, nzt+1
553	u(k,nyn+1:nyng,i) = force%u_north(0,k,i) + ddt_lsf * t_ref * &
554	( force%u_north(1,k,i) - force%u_north(0,k,i) ) * &
555	MERGE( 1.0_wp, 0.0_wp, &
556	BTEST( wall_flags_0(k,nyn+1:nyng,i), 1 ) )
557
558	ENDDO
559	ENDDO
560
561	IF ( .NOT. neutral ) THEN
562	DO i = nxl, nxr
563	DO k = nzb+1, nzt+1
564	pt(k,nyn+1:nyng,i) = force%pt_north(0,k,i) + ddt_lsf * &
565	t_ref * &
566	( force%pt_north(1,k,i) - force%pt_north(0,k,i) )
567
568	ENDDO
569	ENDDO
570	ENDIF
571
572	IF ( humidity ) THEN
573	DO i = nxl, nxr
574	DO k = nzb+1, nzt+1
575	q(k,nyn+1:nyng,i) = force%q_north(0,k,i) + ddt_lsf * &
576	t_ref * &
577	( force%q_north(1,k,i) - force%q_north(0,k,i) )
578
579	ENDDO
580	ENDDO
581	ENDIF
582
583	ENDIF
584	!
585	!-- Top boundary.
586	!-- Please note, only map Inifor data on model top in case the numeric is
587	!-- identical to the Inifor grid. At the top boundary an extrapolation is
588	!-- not possible.
589	DO i = nxlu, nxr
590	DO j = nys, nyn
591	u(nzt+1,j,i) = force%u_top(0,j,i) + ddt_lsf * t_ref * &
592	( force%u_top(1,j,i) - force%u_top(0,j,i) ) * &
593	MERGE( 1.0_wp, 0.0_wp, &
594	BTEST( wall_flags_0(nzt+1,j,i), 1 ) )
595	ENDDO
596	ENDDO
597
598	DO i = nxl, nxr
599	DO j = nysv, nyn
600	v(nzt+1,j,i) = force%v_top(0,j,i) + ddt_lsf * t_ref * &
601	( force%v_top(1,j,i) - force%v_top(0,j,i) ) * &
602	MERGE( 1.0_wp, 0.0_wp, &
603	BTEST( wall_flags_0(nzt+1,j,i), 2 ) )
604	ENDDO
605	ENDDO
606
607	DO i = nxl, nxr
608	DO j = nys, nyn
609	w(nzt:nzt+1,j,i) = force%w_top(0,j,i) + ddt_lsf * t_ref * &
610	( force%w_top(1,j,i) - force%w_top(0,j,i) ) * &
611	MERGE( 1.0_wp, 0.0_wp, &
612	BTEST( wall_flags_0(nzt:nzt+1,j,i), 3 ) )
613	ENDDO
614	ENDDO
615
616
617	IF ( .NOT. neutral ) THEN
618	DO i = nxl, nxr
619	DO j = nys, nyn
620	pt(nzt+1,j,i) = force%pt_top(0,j,i) + ddt_lsf * t_ref * &
621	( force%pt_top(1,j,i) - force%pt_top(0,j,i) )
622	ENDDO
623	ENDDO
624	ENDIF
625
626	IF ( humidity ) THEN
627	DO i = nxl, nxr
628	DO j = nys, nyn
629	q(nzt+1,j,i) = force%q_top(0,j,i) + ddt_lsf * t_ref * &
630	( force%q_top(1,j,i) - force%q_top(0,j,i) )
631	ENDDO
632	ENDDO
633	ENDIF
634	!
635	!-- At the edges( left-south, left-north, right-south and right-north) set
636	!-- data on ghost points.
637	IF ( force_bound_l .AND. force_bound_s ) THEN
638	DO i = 1, nbgp
639	u(:,nys-i,nxlg:nxl) = u(:,nys,nxlg:nxl)
640	w(:,nys-i,nxlg:nxl-1) = w(:,nys,nxlg:nxl-1)
641	IF ( .NOT. neutral ) pt(:,nys-i,nxlg:nxl-1) = pt(:,nys,nxlg:nxl-1)
642	IF ( humidity ) q(:,nys-i,nxlg:nxl-1) = q(:,nys,nxlg:nxl-1)
643	ENDDO
644	DO i = 1, nbgp+1
645	v(:,nysv-i,nxlg:nxl-1) = v(:,nysv,nxlg:nxl-1)
646	ENDDO
647	ENDIF
648	IF ( force_bound_l .AND. force_bound_n ) THEN
649	DO i = 1, nbgp
650	u(:,nyn+i,nxlg:nxl) = u(:,nyn,nxlg:nxl)
651	v(:,nyn+i,nxlg:nxl-1) = v(:,nyn,nxlg:nxl-1)
652	w(:,nyn+i,nxlg:nxl-1) = w(:,nyn,nxlg:nxl-1)
653	IF ( .NOT. neutral ) pt(:,nyn+i,nxlg:nxl-1) = pt(:,nyn,nxlg:nxl-1)
654	IF ( humidity ) q(:,nyn+i,nxlg:nxl-1) = q(:,nyn,nxlg:nxl-1)
655	ENDDO
656	ENDIF
657	IF ( force_bound_r .AND. force_bound_s ) THEN
658	DO i = 1, nbgp
659	u(:,nys-i,nxr+1:nxrg) = u(:,nys,nxr+1:nxrg)
660	w(:,nys-i,nxr+1:nxrg) = w(:,nys,nxr+1:nxrg)
661	IF ( .NOT. neutral ) pt(:,nys-i,nxr+1:nxrg) = pt(:,nys,nxr+1:nxrg)
662	IF ( humidity ) q(:,nys-i,nxr+1:nxrg) = q(:,nys,nxr+1:nxrg)
663	ENDDO
664	DO i = 1, nbgp+1
665	v(:,nysv-i,nxr+1:nxrg) = v(:,nysv,nxr+1:nxrg)
666	ENDDO
667	ENDIF
668	IF ( force_bound_r .AND. force_bound_n ) THEN
669	DO i = 1, nbgp
670	u(:,nyn+i,nxr+1:nxrg) = u(:,nyn,nxr+1:nxrg)
671	v(:,nyn+i,nxr+1:nxrg) = v(:,nyn,nxr+1:nxrg)
672	w(:,nyn+i,nxr+1:nxrg) = w(:,nyn,nxr+1:nxrg)
673	IF ( .NOT. neutral ) pt(:,nyn+i,nxr+1:nxrg) = pt(:,nyn,nxr+1:nxrg)
674	IF ( humidity ) q(:,nyn+i,nxr+1:nxrg) = q(:,nyn,nxr+1:nxrg)
675	ENDDO
676	ENDIF
677	!
678	!-- Moreover, set Neumann boundary condition for subgrid-scale TKE and
679	!-- passive scalar
680	IF ( .NOT. constant_diffusion ) THEN
681	IF ( force_bound_l ) e(:,:,nxl-1) = e(:,:,nxl)
682	IF ( force_bound_r ) e(:,:,nxr+1) = e(:,:,nxr)
683	IF ( force_bound_s ) e(:,nys-1,:) = e(:,nys,:)
684	IF ( force_bound_n ) e(:,nyn+1,:) = e(:,nyn,:)
685	e(nzt+1,:,:) = e(nzt,:,:)
686	ENDIF
687	IF ( passive_scalar ) THEN
688	IF ( force_bound_l ) s(:,:,nxl-1) = s(:,:,nxl)
689	IF ( force_bound_r ) s(:,:,nxr+1) = s(:,:,nxr)
690	IF ( force_bound_s ) s(:,nys-1,:) = s(:,nys,:)
691	IF ( force_bound_n ) s(:,nyn+1,:) = s(:,nyn,:)
692	ENDIF
693
694
695	CALL exchange_horiz( u, nbgp )
696	CALL exchange_horiz( v, nbgp )
697	CALL exchange_horiz( w, nbgp )
698	IF ( .NOT. neutral ) CALL exchange_horiz( pt, nbgp )
699	IF ( humidity ) CALL exchange_horiz( q, nbgp )
700
701	!
702	!-- Set surface pressure. Please note, time-dependent surface
703	!-- pressure would require changes in anelastic approximation and
704	!-- treatment of fluxes.
705	!-- For the moment, comment this out!
706	! surface_pressure = force%surface_pressure(force%tind) + &
707	! ddt_lsf * t_ref * &
708	! ( force%surface_pressure(force%tind_p) &
709	! - force%surface_pressure(force%tind) )
710
711	END SUBROUTINE forcing_bc
712
713	!------------------------------------------------------------------------------!
714	! Description:
715	! ------------
716	!> @todo Missing subroutine description.
717	!------------------------------------------------------------------------------!
718	SUBROUTINE lsf_nudging_check_parameters
719
720	IMPLICIT NONE
721	!
722	!-- Check nudging and large scale forcing from external file
723	IF ( nudging .AND. ( .NOT. large_scale_forcing ) ) THEN
724	message_string = 'Nudging requires large_scale_forcing = .T.. &'// &
725	'Surface fluxes and geostrophic wind should be &'// &
726	'prescribed in file LSF_DATA'
727	CALL message( 'check_parameters', 'PA0374', 1, 2, 0, 6, 0 )
728	ENDIF
729
730	IF ( large_scale_forcing .AND. ( bc_lr /= 'cyclic' .OR. &
731	bc_ns /= 'cyclic' ) ) THEN
732	message_string = 'Non-cyclic lateral boundaries do not allow for &'//&
733	'the usage of large scale forcing from external file.'
734	CALL message( 'check_parameters', 'PA0375', 1, 2, 0, 6, 0 )
735	ENDIF
736
737	IF ( large_scale_forcing .AND. ( .NOT. humidity ) ) THEN
738	message_string = 'The usage of large scale forcing from external &'//&
739	'file LSF_DATA requires humidity = .T..'
740	CALL message( 'check_parameters', 'PA0376', 1, 2, 0, 6, 0 )
741	ENDIF
742
743	IF ( large_scale_forcing .AND. passive_scalar ) THEN
744	message_string = 'The usage of large scale forcing from external &'//&
745	'file LSF_DATA is not implemented for passive scalars'
746	CALL message( 'check_parameters', 'PA0440', 1, 2, 0, 6, 0 )
747	ENDIF
748
749	IF ( large_scale_forcing .AND. topography /= 'flat' &
750	.AND. .NOT. lsf_exception ) THEN
751	message_string = 'The usage of large scale forcing from external &'//&
752	'file LSF_DATA is not implemented for non-flat topography'
753	CALL message( 'check_parameters', 'PA0377', 1, 2, 0, 6, 0 )
754	ENDIF
755
756	IF ( large_scale_forcing .AND. ocean ) THEN
757	message_string = 'The usage of large scale forcing from external &'//&
758	'file LSF_DATA is not implemented for ocean runs'
759	CALL message( 'check_parameters', 'PA0378', 1, 2, 0, 6, 0 )
760	ENDIF
761
762	END SUBROUTINE lsf_nudging_check_parameters
763
764	!------------------------------------------------------------------------------!
765	! Description:
766	! ------------
767	!> Check data output of profiles for land surface model
768	!------------------------------------------------------------------------------!
769	SUBROUTINE lsf_nudging_check_data_output_pr( variable, var_count, unit, &
770	dopr_unit )
771
772	USE profil_parameter
773
774	IMPLICIT NONE
775
776	CHARACTER (LEN=*) :: unit !<
777	CHARACTER (LEN=*) :: variable !<
778	CHARACTER (LEN=*) :: dopr_unit !< local value of dopr_unit
779
780	INTEGER(iwp) :: user_pr_index !<
781	INTEGER(iwp) :: var_count !<
782
783	SELECT CASE ( TRIM( variable ) )
784
785
786	CASE ( 'td_lsa_lpt' )
787	IF ( .NOT. large_scale_forcing ) THEN
788	message_string = 'data_output_pr = ' // &
789	TRIM( data_output_pr(var_count) ) // &
790	' is not implemented for ' // &
791	'large_scale_forcing = .FALSE.'
792	CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393', &
793	1, 2, 0, 6, 0 )
794	ELSE
795	dopr_index(var_count) = 81
796	dopr_unit = 'K/s'
797	unit = 'K/s'
798	hom(:,2,81,:) = SPREAD( zu, 2, statistic_regions+1 )
799	ENDIF
800
801	CASE ( 'td_lsa_q' )
802	IF ( .NOT. large_scale_forcing ) THEN
803	message_string = 'data_output_pr = ' // &
804	TRIM( data_output_pr(var_count) ) // &
805	' is not implemented for ' // &
806	'large_scale_forcing = .FALSE.'
807	CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393', &
808	1, 2, 0, 6, 0 )
809	ELSE
810	dopr_index(var_count) = 82
811	dopr_unit = 'kg/kgs'
812	unit = 'kg/kgs'
813	hom(:,2,82,:) = SPREAD( zu, 2, statistic_regions+1 )
814	ENDIF
815	CASE ( 'td_sub_lpt' )
816	IF ( .NOT. large_scale_forcing ) THEN
817	message_string = 'data_output_pr = ' // &
818	TRIM( data_output_pr(var_count) ) // &
819	' is not implemented for ' // &
820	'large_scale_forcing = .FALSE.'
821	CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393', &
822	1, 2, 0, 6, 0 )
823	ELSE
824	dopr_index(var_count) = 83
825	dopr_unit = 'K/s'
826	unit = 'K/s'
827	hom(:,2,83,:) = SPREAD( zu, 2, statistic_regions+1 )
828	ENDIF
829
830	CASE ( 'td_sub_q' )
831	IF ( .NOT. large_scale_forcing ) THEN
832	message_string = 'data_output_pr = ' // &
833	TRIM( data_output_pr(var_count) ) // &
834	' is not implemented for ' // &
835	'large_scale_forcing = .FALSE.'
836	CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393', &
837	1, 2, 0, 6, 0 )
838	ELSE
839	dopr_index(var_count) = 84
840	dopr_unit = 'kg/kgs'
841	unit = 'kg/kgs'
842	hom(:,2,84,:) = SPREAD( zu, 2, statistic_regions+1 )
843	ENDIF
844
845	CASE ( 'td_nud_lpt' )
846	IF ( .NOT. nudging ) THEN
847	message_string = 'data_output_pr = ' // &
848	TRIM( data_output_pr(var_count) ) // &
849	' is not implemented for ' // &
850	'large_scale_forcing = .FALSE.'
851	CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393', &
852	1, 2, 0, 6, 0 )
853	ELSE
854	dopr_index(var_count) = 85
855	dopr_unit = 'K/s'
856	unit = 'K/s'
857	hom(:,2,85,:) = SPREAD( zu, 2, statistic_regions+1 )
858	ENDIF
859
860	CASE ( 'td_nud_q' )
861	IF ( .NOT. nudging ) THEN
862	message_string = 'data_output_pr = ' // &
863	TRIM( data_output_pr(var_count) ) // &
864	' is not implemented for ' // &
865	'large_scale_forcing = .FALSE.'
866	CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393', &
867	1, 2, 0, 6, 0 )
868	ELSE
869	dopr_index(var_count) = 86
870	dopr_unit = 'kg/kgs'
871	unit = 'kg/kgs'
872	hom(:,2,86,:) = SPREAD( zu, 2, statistic_regions+1 )
873	ENDIF
874
875	CASE ( 'td_nud_u' )
876	IF ( .NOT. nudging ) THEN
877	message_string = 'data_output_pr = ' // &
878	TRIM( data_output_pr(var_count) ) // &
879	' is not implemented for ' // &
880	'large_scale_forcing = .FALSE.'
881	CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393', &
882	1, 2, 0, 6, 0 )
883	ELSE
884	dopr_index(var_count) = 87
885	dopr_unit = 'm/s2'
886	unit = 'm/s2'
887	hom(:,2,87,:) = SPREAD( zu, 2, statistic_regions+1 )
888	ENDIF
889
890	CASE ( 'td_nud_v' )
891	IF ( .NOT. nudging ) THEN
892	message_string = 'data_output_pr = ' // &
893	TRIM( data_output_pr(var_count) ) // &
894	' is not implemented for ' // &
895	'large_scale_forcing = .FALSE.'
896	CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393', &
897	1, 2, 0, 6, 0 )
898	ELSE
899	dopr_index(var_count) = 88
900	dopr_unit = 'm/s2'
901	unit = 'm/s2'
902	hom(:,2,88,:) = SPREAD( zu, 2, statistic_regions+1 )
903	ENDIF
904
905
906	CASE DEFAULT
907	unit = 'illegal'
908
909	END SELECT
910
911	END SUBROUTINE lsf_nudging_check_data_output_pr
912
913	!------------------------------------------------------------------------------!
914	! Description:
915	! ------------
916	!> @todo Missing subroutine description.
917	!------------------------------------------------------------------------------!
918	SUBROUTINE lsf_nudging_header ( io )
919
920	IMPLICIT NONE
921
922	INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file
923
924	WRITE ( io, 1 )
925	IF ( large_scale_forcing ) THEN
926	WRITE ( io, 3 )
927	WRITE ( io, 4 )
928
929	IF ( large_scale_subsidence ) THEN
930	IF ( .NOT. use_subsidence_tendencies ) THEN
931	WRITE ( io, 5 )
932	ELSE
933	WRITE ( io, 6 )
934	ENDIF
935	ENDIF
936
937	IF ( bc_pt_b == 'dirichlet' ) THEN
938	WRITE ( io, 12 )
939	ELSEIF ( bc_pt_b == 'neumann' ) THEN
940	WRITE ( io, 13 )
941	ENDIF
942
943	IF ( bc_q_b == 'dirichlet' ) THEN
944	WRITE ( io, 14 )
945	ELSEIF ( bc_q_b == 'neumann' ) THEN
946	WRITE ( io, 15 )
947	ENDIF
948
949	WRITE ( io, 7 )
950	IF ( nudging ) THEN
951	WRITE ( io, 10 )
952	ENDIF
953	ELSE
954	WRITE ( io, 2 )
955	WRITE ( io, 11 )
956	ENDIF
957	IF ( large_scale_subsidence ) THEN
958	WRITE ( io, 8 )
959	WRITE ( io, 9 )
960	ENDIF
961
962
963	1 FORMAT (//' Large scale forcing and nudging:'/ &
964	' -------------------------------'/)
965	2 FORMAT (' --> No large scale forcing from external is used (default) ')
966	3 FORMAT (' --> Large scale forcing from external file LSF_DATA is used: ')
967	4 FORMAT (' - large scale advection tendencies ')
968	5 FORMAT (' - large scale subsidence velocity w_subs ')
969	6 FORMAT (' - large scale subsidence tendencies ')
970	7 FORMAT (' - and geostrophic wind components ug and vg')
971	8 FORMAT (' --> Large-scale vertical motion is used in the ', &
972	'prognostic equation(s) for')
973	9 FORMAT (' the scalar(s) only')
974	10 FORMAT (' --> Nudging is used')
975	11 FORMAT (' --> No nudging is used (default) ')
976	12 FORMAT (' - prescribed surface values for temperature')
977	13 FORMAT (' - prescribed surface fluxes for temperature')
978	14 FORMAT (' - prescribed surface values for humidity')
979	15 FORMAT (' - prescribed surface fluxes for humidity')
980
981	END SUBROUTINE lsf_nudging_header
982
983	!------------------------------------------------------------------------------!
984	! Description:
985	! ------------
986	!> @todo Missing subroutine description.
987	!------------------------------------------------------------------------------!
988	SUBROUTINE lsf_init
989
990	USE control_parameters, &
991	ONLY: bc_lr_cyc, bc_ns_cyc
992
993	USE netcdf_data_input_mod, &
994	ONLY: netcdf_data_input_lsf
995
996	IMPLICIT NONE
997
998	CHARACTER(100) :: chmess !<
999	CHARACTER(1) :: hash !<
1000
1001	INTEGER(iwp) :: ierrn !<
1002	INTEGER(iwp) :: finput = 90 !<
1003	INTEGER(iwp) :: k !<
1004	INTEGER(iwp) :: nt !<
1005	INTEGER(iwp) :: t !< running index for time levels
1006
1007	REAL(wp) :: fac !<
1008	REAL(wp) :: highheight !<
1009	REAL(wp) :: highug_vert !<
1010	REAL(wp) :: highvg_vert !<
1011	REAL(wp) :: highwsubs_vert !<
1012	REAL(wp) :: lowheight !<
1013	REAL(wp) :: lowug_vert !<
1014	REAL(wp) :: lowvg_vert !<
1015	REAL(wp) :: lowwsubs_vert !<
1016	REAL(wp) :: high_td_lsa_lpt !<
1017	REAL(wp) :: low_td_lsa_lpt !<
1018	REAL(wp) :: high_td_lsa_q !<
1019	REAL(wp) :: low_td_lsa_q !<
1020	REAL(wp) :: high_td_sub_lpt !<
1021	REAL(wp) :: low_td_sub_lpt !<
1022	REAL(wp) :: high_td_sub_q !<
1023	REAL(wp) :: low_td_sub_q !<
1024	REAL(wp) :: r_dummy !<
1025
1026	IF ( forcing ) THEN
1027	!
1028	!-- Allocate arrays for geostrophic wind components. Arrays will
1029	!-- incorporate 2 time levels in order to interpolate in between. Please
1030	!-- note, forcing using geostrophic wind components is only required in
1031	!-- case of cyclic boundary conditions.
1032	IF ( bc_lr_cyc .AND. bc_ns_cyc ) THEN
1033	ALLOCATE( force%ug(0:1,nzb:nzt+1) )
1034	ALLOCATE( force%vg(0:1,nzb:nzt+1) )
1035	ENDIF
1036	!
1037	!-- Allocate arrays for reading boundary values. Arrays will incorporate 2
1038	!-- time levels in order to interpolate in between.
1039	IF ( force_bound_l ) THEN
1040	ALLOCATE( force%u_left(0:1,nzb+1:nzt+1,nys:nyn) )
1041	ALLOCATE( force%v_left(0:1,nzb+1:nzt+1,nysv:nyn) )
1042	ALLOCATE( force%w_left(0:1,nzb+1:nzt,nys:nyn) )
1043	IF ( humidity ) ALLOCATE( force%q_left(0:1,nzb+1:nzt+1,nys:nyn) )
1044	IF ( .NOT. neutral ) ALLOCATE( force%pt_left(0:1,nzb+1:nzt+1,nys:nyn) )
1045	ENDIF
1046	IF ( force_bound_r ) THEN
1047	ALLOCATE( force%u_right(0:1,nzb+1:nzt+1,nys:nyn) )
1048	ALLOCATE( force%v_right(0:1,nzb+1:nzt+1,nysv:nyn) )
1049	ALLOCATE( force%w_right(0:1,nzb+1:nzt,nys:nyn) )
1050	IF ( humidity ) ALLOCATE( force%q_right(0:1,nzb+1:nzt+1,nys:nyn) )
1051	IF ( .NOT. neutral ) ALLOCATE( force%pt_right(0:1,nzb+1:nzt+1,nys:nyn) )
1052	ENDIF
1053	IF ( force_bound_n ) THEN
1054	ALLOCATE( force%u_north(0:1,nzb+1:nzt+1,nxlu:nxr) )
1055	ALLOCATE( force%v_north(0:1,nzb+1:nzt+1,nxl:nxr) )
1056	ALLOCATE( force%w_north(0:1,nzb+1:nzt,nxl:nxr) )
1057	IF ( humidity ) ALLOCATE( force%q_north(0:1,nzb+1:nzt+1,nxl:nxr) )
1058	IF ( .NOT. neutral ) ALLOCATE( force%pt_north(0:1,nzb+1:nzt+1,nxl:nxr) )
1059	ENDIF
1060	IF ( force_bound_s ) THEN
1061	ALLOCATE( force%u_south(0:1,nzb+1:nzt+1,nxlu:nxr) )
1062	ALLOCATE( force%v_south(0:1,nzb+1:nzt+1,nxl:nxr) )
1063	ALLOCATE( force%w_south(0:1,nzb+1:nzt,nxl:nxr) )
1064	IF ( humidity ) ALLOCATE( force%q_south(0:1,nzb+1:nzt+1,nxl:nxr) )
1065	IF ( .NOT. neutral ) ALLOCATE( force%pt_south(0:1,nzb+1:nzt+1,nxl:nxr) )
1066	ENDIF
1067
1068	ALLOCATE( force%u_top(0:1,nys:nyn,nxlu:nxr) )
1069	ALLOCATE( force%v_top(0:1,nysv:nyn,nxl:nxr) )
1070	ALLOCATE( force%w_top(0:1,nys:nyn,nxl:nxr) )
1071	IF ( humidity ) ALLOCATE( force%q_top(0:1,nys:nyn,nxl:nxr) )
1072	IF ( .NOT. neutral ) ALLOCATE( force%pt_top(0:1,nys:nyn,nxl:nxr) )
1073
1074	!
1075	!-- Initial call of input. Time array, initial 3D data of u, v, w,
1076	!-- potential temperature, as well as mixing ratio, will be read.
1077	!-- Moreover, data at lateral and top boundary will be read.
1078	CALL netcdf_data_input_lsf
1079	!
1080	!-- Please note, at the moment INIFOR assumes only an equidistant vertical
1081	!-- grid. In case of vertical grid stretching, input of inital 3D data
1082	!-- need to be inter- and/or extrapolated.
1083	!-- Therefore, check if zw grid on file is identical to numeric zw grid.
1084	IF ( ANY( zu(1:nzt+1) /= force%zu_atmos(1:force%nzu) ) ) THEN
1085	!
1086	!-- Also data at the boundaries need to be inter/extrapolated at both
1087	!-- time levels
1088	DO t = 0, 1
1089	IF ( force_bound_l ) THEN
1090	CALL netcdf_data_input_interpolate( force%u_left(t,:,:), &
1091	zu(1:nzt+1), &
1092	force%zu_atmos )
1093	CALL netcdf_data_input_interpolate( force%v_left(t,:,:), &
1094	zu(1:nzt+1), &
1095	force%zu_atmos )
1096	CALL netcdf_data_input_interpolate( force%w_left(t,:,:), &
1097	zw(1:nzt+1), &
1098	force%zw_atmos )
1099	IF ( .NOT. neutral ) &
1100	CALL netcdf_data_input_interpolate( force%pt_left(t,:,:),&
1101	zu(1:nzt+1), &
1102	force%zu_atmos )
1103	IF ( humidity ) &
1104	CALL netcdf_data_input_interpolate( force%q_left(t,:,:), &
1105	zu(1:nzt+1), &
1106	force%zu_atmos )
1107	ENDIF
1108	IF ( force_bound_r ) THEN
1109	CALL netcdf_data_input_interpolate( force%u_right(t,:,:), &
1110	zu(1:nzt+1), &
1111	force%zu_atmos )
1112	CALL netcdf_data_input_interpolate( force%v_right(t,:,:), &
1113	zu(1:nzt+1), &
1114	force%zu_atmos )
1115	CALL netcdf_data_input_interpolate( force%w_right(t,:,:), &
1116	zw(1:nzt+1), &
1117	force%zw_atmos )
1118	IF ( .NOT. neutral ) &
1119	CALL netcdf_data_input_interpolate( force%pt_right(t,:,:),&
1120	zu(1:nzt+1), &
1121	force%zu_atmos )
1122	IF ( humidity ) &
1123	CALL netcdf_data_input_interpolate( force%q_right(t,:,:),&
1124	zu(1:nzt+1), &
1125	force%zu_atmos )
1126	ENDIF
1127	IF ( force_bound_n ) THEN
1128	CALL netcdf_data_input_interpolate( force%u_north(t,:,:), &
1129	zu(1:nzt+1), &
1130	force%zu_atmos )
1131	CALL netcdf_data_input_interpolate( force%v_north(t,:,:), &
1132	zu(1:nzt+1), &
1133	force%zu_atmos )
1134	CALL netcdf_data_input_interpolate( force%w_north(t,:,:), &
1135	zw(1:nzt+1), &
1136	force%zw_atmos )
1137	IF ( .NOT. neutral ) &
1138	CALL netcdf_data_input_interpolate( force%pt_north(t,:,:),&
1139	zu(1:nzt+1), &
1140	force%zu_atmos )
1141	IF ( humidity ) &
1142	CALL netcdf_data_input_interpolate( force%q_north(t,:,:),&
1143	zu(1:nzt+1), &
1144	force%zu_atmos )
1145	ENDIF
1146	IF ( force_bound_s ) THEN
1147	CALL netcdf_data_input_interpolate( force%u_south(t,:,:), &
1148	zu(1:nzt+1), &
1149	force%zu_atmos )
1150	CALL netcdf_data_input_interpolate( force%v_south(t,:,:), &
1151	zu(1:nzt+1), &
1152	force%zu_atmos )
1153	CALL netcdf_data_input_interpolate( force%w_south(t,:,:), &
1154	zw(1:nzt+1), &
1155	force%zw_atmos )
1156	IF ( .NOT. neutral ) &
1157	CALL netcdf_data_input_interpolate( force%pt_south(t,:,:),&
1158	zu(1:nzt+1), &
1159	force%zu_atmos )
1160	IF ( humidity ) &
1161	CALL netcdf_data_input_interpolate( force%q_south(t,:,:),&
1162	zu(1:nzt+1), &
1163	force%zu_atmos )
1164	ENDIF
1165	ENDDO
1166	ENDIF
1167
1168	!
1169	!-- Exchange ghost points
1170	CALL exchange_horiz( u, nbgp )
1171	CALL exchange_horiz( v, nbgp )
1172	CALL exchange_horiz( w, nbgp )
1173	IF ( .NOT. neutral ) CALL exchange_horiz( pt, nbgp )
1174	IF ( humidity ) CALL exchange_horiz( q, nbgp )
1175	!
1176	!-- At lateral boundaries, set also initial boundary conditions
1177	IF ( force_bound_l ) THEN
1178	u(:,:,nxl) = u(:,:,nxlu)
1179	v(:,:,nxl-1) = v(:,:,nxl)
1180	w(:,:,nxl-1) = w(:,:,nxl)
1181	IF ( .NOT. neutral ) pt(:,:,nxl-1) = pt(:,:,nxl)
1182	IF ( humidity ) q(:,:,nxl-1) = q(:,:,nxl)
1183	ENDIF
1184	IF ( force_bound_r ) THEN
1185	u(:,:,nxr+1) = u(:,:,nxr)
1186	v(:,:,nxr+1) = v(:,:,nxr)
1187	w(:,:,nxr+1) = w(:,:,nxr)
1188	IF ( .NOT. neutral ) pt(:,:,nxr+1) = pt(:,:,nxr)
1189	IF ( humidity ) q(:,:,nxr+1) = q(:,:,nxr)
1190	ENDIF
1191	IF ( force_bound_s ) THEN
1192	u(:,nys-1,:) = u(:,nys,:)
1193	v(:,nys,:) = v(:,nysv,:)
1194	w(:,nys-1,:) = w(:,nys,:)
1195	IF ( .NOT. neutral ) pt(:,nys-1,:) = pt(:,nys,:)
1196	IF ( humidity ) q(:,nys-1,:) = q(:,nys,:)
1197	ENDIF
1198	IF ( force_bound_n ) THEN
1199	u(:,nyn+1,:) = u(:,nyn,:)
1200	v(:,nyn+1,:) = v(:,nyn,:)
1201	w(:,nyn+1,:) = w(:,nyn,:)
1202	IF ( .NOT. neutral ) pt(:,nyn+1,:) = pt(:,nyn,:)
1203	IF ( humidity ) q(:,nyn+1,:) = q(:,nyn,:)
1204	ENDIF
1205
1206	!
1207	!-- After 3D data is initialized, ensure mass conservation
1208	CALL forcing_bc_mass_conservation
1209	!
1210	!-- Initialize surface pressure. Please note, time-dependent surface
1211	!-- pressure would require changes in anelastic approximation and
1212	!-- treatment of fluxes.
1213	!-- For the moment, comment this out!
1214	! surface_pressure = force%surface_pressure(0)
1215
1216	ELSE
1217
1218	ALLOCATE( p_surf(0:nlsf), pt_surf(0:nlsf), q_surf(0:nlsf), &
1219	qsws_surf(0:nlsf), shf_surf(0:nlsf), &
1220	td_lsa_lpt(nzb:nzt+1,0:nlsf), td_lsa_q(nzb:nzt+1,0:nlsf), &
1221	td_sub_lpt(nzb:nzt+1,0:nlsf), td_sub_q(nzb:nzt+1,0:nlsf), &
1222	time_vert(0:nlsf), time_surf(0:nlsf), &
1223	ug_vert(nzb:nzt+1,0:nlsf), vg_vert(nzb:nzt+1,0:nlsf), &
1224	wsubs_vert(nzb:nzt+1,0:nlsf) )
1225
1226	p_surf = 0.0_wp; pt_surf = 0.0_wp; q_surf = 0.0_wp; qsws_surf = 0.0_wp
1227	shf_surf = 0.0_wp; time_vert = 0.0_wp; td_lsa_lpt = 0.0_wp
1228	td_lsa_q = 0.0_wp; td_sub_lpt = 0.0_wp; td_sub_q = 0.0_wp
1229	time_surf = 0.0_wp; ug_vert = 0.0_wp; vg_vert = 0.0_wp
1230	wsubs_vert = 0.0_wp
1231
1232	!
1233	!-- Array for storing large scale forcing and nudging tendencies at each
1234	!-- timestep for data output
1235	ALLOCATE( sums_ls_l(nzb:nzt+1,0:7) )
1236	sums_ls_l = 0.0_wp
1237
1238	ngp_sums_ls = (nz+2)*6
1239
1240	OPEN ( finput, FILE='LSF_DATA', STATUS='OLD', &
1241	FORM='FORMATTED', IOSTAT=ierrn )
1242
1243	IF ( ierrn /= 0 ) THEN
1244	message_string = 'file LSF_DATA does not exist'
1245	CALL message( 'ls_forcing', 'PA0368', 1, 2, 0, 6, 0 )
1246	ENDIF
1247
1248	ierrn = 0
1249	!
1250	!-- First three lines of LSF_DATA contain header
1251	READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
1252	READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
1253	READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
1254
1255	IF ( ierrn /= 0 ) THEN
1256	message_string = 'errors in file LSF_DATA'
1257	CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
1258	ENDIF
1259
1260	!
1261	!-- Surface values are read in
1262	nt = 0
1263	ierrn = 0
1264
1265	DO WHILE ( time_surf(nt) < end_time )
1266	nt = nt + 1
1267	READ ( finput, *, IOSTAT = ierrn ) time_surf(nt), shf_surf(nt), &
1268	qsws_surf(nt), pt_surf(nt), &
1269	q_surf(nt), p_surf(nt)
1270
1271	IF ( ierrn /= 0 ) THEN
1272	WRITE ( message_string, * ) 'No time dependent surface variables ' //&
1273	'in&LSF_DATA for end of run found'
1274
1275	CALL message( 'ls_forcing', 'PA0363', 1, 2, 0, 6, 0 )
1276	ENDIF
1277	ENDDO
1278
1279	IF ( time_surf(1) > end_time ) THEN
1280	WRITE ( message_string, * ) 'Time dependent surface variables in ' // &
1281	'&LSF_DATA set in after end of ' , &
1282	'simulation - lsf_surf is set to FALSE'
1283	CALL message( 'ls_forcing', 'PA0371', 0, 0, 0, 6, 0 )
1284	lsf_surf = .FALSE.
1285	ENDIF
1286
1287	!
1288	!-- Go to the end of the list with surface variables
1289	DO WHILE ( ierrn == 0 )
1290	READ ( finput, *, IOSTAT = ierrn ) r_dummy
1291	ENDDO
1292
1293	!
1294	!-- Profiles of ug, vg and w_subs are read in (large scale forcing)
1295
1296	nt = 0
1297	DO WHILE ( time_vert(nt) < end_time )
1298	nt = nt + 1
1299	hash = "#"
1300	ierrn = 1 ! not zero
1301	!
1302	!-- Search for the next line consisting of "# time",
1303	!-- from there onwards the profiles will be read
1304	DO WHILE ( .NOT. ( hash == "#" .AND. ierrn == 0 ) )
1305	READ ( finput, *, IOSTAT=ierrn ) hash, time_vert(nt)
1306	IF ( ierrn < 0 ) THEN
1307	WRITE( message_string, * ) 'No time dependent vertical profiles',&
1308	' in&LSF_DATA for end of run found'
1309	CALL message( 'ls_forcing', 'PA0372', 1, 2, 0, 6, 0 )
1310	ENDIF
1311	ENDDO
1312
1313	IF ( nt == 1 .AND. time_vert(nt) > end_time ) EXIT
1314
1315	READ ( finput, *, IOSTAT=ierrn ) lowheight, lowug_vert, lowvg_vert,&
1316	lowwsubs_vert, low_td_lsa_lpt, &
1317	low_td_lsa_q, low_td_sub_lpt, &
1318	low_td_sub_q
1319	IF ( ierrn /= 0 ) THEN
1320	message_string = 'errors in file LSF_DATA'
1321	CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
1322	ENDIF
1323
1324	READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert, &
1325	highvg_vert, highwsubs_vert, &
1326	high_td_lsa_lpt, high_td_lsa_q, &
1327	high_td_sub_lpt, high_td_sub_q
1328
1329	IF ( ierrn /= 0 ) THEN
1330	message_string = 'errors in file LSF_DATA'
1331	CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
1332	ENDIF
1333
1334
1335	DO k = nzb, nzt+1
1336	IF ( highheight < zu(k) ) THEN
1337	lowheight = highheight
1338	lowug_vert = highug_vert
1339	lowvg_vert = highvg_vert
1340	lowwsubs_vert = highwsubs_vert
1341	low_td_lsa_lpt = high_td_lsa_lpt
1342	low_td_lsa_q = high_td_lsa_q
1343	low_td_sub_lpt = high_td_sub_lpt
1344	low_td_sub_q = high_td_sub_q
1345
1346	ierrn = 0
1347	READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert, &
1348	highvg_vert, highwsubs_vert,&
1349	high_td_lsa_lpt, &
1350	high_td_lsa_q, &
1351	high_td_sub_lpt, high_td_sub_q
1352
1353	IF ( ierrn /= 0 ) THEN
1354	WRITE( message_string, * ) 'zu(',k,') = ', zu(k), 'm ', &
1355	'is higher than the maximum height in LSF_DATA which ',&
1356	'is ', lowheight, 'm. Interpolation on PALM ', &
1357	'grid is not possible.'
1358	CALL message( 'ls_forcing', 'PA0395', 1, 2, 0, 6, 0 )
1359	ENDIF
1360
1361	ENDIF
1362
1363	!
1364	!-- Interpolation of prescribed profiles in space
1365	fac = (highheight-zu(k))/(highheight - lowheight)
1366
1367	ug_vert(k,nt) = fac * lowug_vert &
1368	+ ( 1.0_wp - fac ) * highug_vert
1369	vg_vert(k,nt) = fac * lowvg_vert &
1370	+ ( 1.0_wp - fac ) * highvg_vert
1371	wsubs_vert(k,nt) = fac * lowwsubs_vert &
1372	+ ( 1.0_wp - fac ) * highwsubs_vert
1373
1374	td_lsa_lpt(k,nt) = fac * low_td_lsa_lpt &
1375	+ ( 1.0_wp - fac ) * high_td_lsa_lpt
1376	td_lsa_q(k,nt) = fac * low_td_lsa_q &
1377	+ ( 1.0_wp - fac ) * high_td_lsa_q
1378	td_sub_lpt(k,nt) = fac * low_td_sub_lpt &
1379	+ ( 1.0_wp - fac ) * high_td_sub_lpt
1380	td_sub_q(k,nt) = fac * low_td_sub_q &
1381	+ ( 1.0_wp - fac ) * high_td_sub_q
1382
1383	ENDDO
1384
1385	ENDDO
1386
1387	!
1388	!-- Large scale vertical velocity has to be zero at the surface
1389	wsubs_vert(nzb,:) = 0.0_wp
1390
1391	IF ( time_vert(1) > end_time ) THEN
1392	WRITE ( message_string, * ) 'Time dependent large scale profile ',&
1393	'forcing from&LSF_DATA sets in after end of ' ,&
1394	'simulation - lsf_vert is set to FALSE'
1395	CALL message( 'ls_forcing', 'PA0373', 0, 0, 0, 6, 0 )
1396	lsf_vert = .FALSE.
1397	ENDIF
1398
1399	CLOSE( finput )
1400
1401	ENDIF
1402
1403	END SUBROUTINE lsf_init
1404
1405	!------------------------------------------------------------------------------!
1406	! Description:
1407	! ------------
1408	!> @todo Missing subroutine description.
1409	!------------------------------------------------------------------------------!
1410	SUBROUTINE ls_forcing_surf ( time )
1411
1412	IMPLICIT NONE
1413
1414	INTEGER(iwp) :: nt !<
1415
1416	REAL(wp) :: dum_surf_flux !<
1417	REAL(wp) :: fac !<
1418	REAL(wp), INTENT(in) :: time !<
1419
1420	!
1421	!-- Interpolation in time of LSF_DATA at the surface
1422	nt = 1
1423	DO WHILE ( time > time_surf(nt) )
1424	nt = nt + 1
1425	ENDDO
1426	IF ( time /= time_surf(nt) ) THEN
1427	nt = nt - 1
1428	ENDIF
1429
1430	fac = ( time -time_surf(nt) ) / ( time_surf(nt+1) - time_surf(nt) )
1431
1432	IF ( ibc_pt_b == 0 ) THEN
1433	!
1434	!-- In case of Dirichlet boundary condition shf must not
1435	!-- be set - it is calculated via MOST in prandtl_fluxes
1436	pt_surface = pt_surf(nt) + fac * ( pt_surf(nt+1) - pt_surf(nt) )
1437
1438	ELSEIF ( ibc_pt_b == 1 ) THEN
1439	!
1440	!-- In case of Neumann boundary condition pt_surface is needed for
1441	!-- calculation of reference density
1442	dum_surf_flux = ( shf_surf(nt) + fac * &
1443	( shf_surf(nt+1) - shf_surf(nt) ) &
1444	) * heatflux_input_conversion(nzb)
1445	!
1446	!-- Save surface sensible heat flux on default, natural and urban surface
1447	!-- type, if required
1448	IF ( surf_def_h(0)%ns >= 1 ) surf_def_h(0)%shf(:) = dum_surf_flux
1449	IF ( surf_lsm_h%ns >= 1 ) surf_lsm_h%shf(:) = dum_surf_flux
1450	IF ( surf_usm_h%ns >= 1 ) surf_usm_h%shf(:) = dum_surf_flux
1451
1452	pt_surface = pt_surf(nt) + fac * ( pt_surf(nt+1) - pt_surf(nt) )
1453
1454	ENDIF
1455
1456	IF ( ibc_q_b == 0 ) THEN
1457	!
1458	!-- In case of Dirichlet boundary condition qsws must not
1459	!-- be set - it is calculated via MOST in prandtl_fluxes
1460	q_surface = q_surf(nt) + fac * ( q_surf(nt+1) - q_surf(nt) )
1461
1462	ELSEIF ( ibc_q_b == 1 ) THEN
1463	dum_surf_flux = ( qsws_surf(nt) + fac * &
1464	( qsws_surf(nt+1) - qsws_surf(nt) ) &
1465	) * waterflux_input_conversion(nzb)
1466	!
1467	!-- Save surface sensible heat flux on default, natural and urban surface
1468	!-- type, if required
1469	IF ( surf_def_h(0)%ns >= 1 ) surf_def_h(0)%qsws(:) = dum_surf_flux
1470	IF ( surf_lsm_h%ns >= 1 ) surf_lsm_h%qsws(:) = dum_surf_flux
1471	IF ( surf_usm_h%ns >= 1 ) surf_usm_h%qsws(:) = dum_surf_flux
1472
1473	ENDIF
1474
1475	surface_pressure = p_surf(nt) + fac * ( p_surf(nt+1) - p_surf(nt) )
1476
1477	END SUBROUTINE ls_forcing_surf
1478
1479
1480
1481
1482	!------------------------------------------------------------------------------!
1483	! Description:
1484	! ------------
1485	!> @todo Missing subroutine description.
1486	!------------------------------------------------------------------------------!
1487	SUBROUTINE ls_forcing_vert ( time )
1488
1489
1490	IMPLICIT NONE
1491
1492	INTEGER(iwp) :: nt !<
1493
1494	REAL(wp) :: fac !<
1495	REAL(wp), INTENT(in) :: time !<
1496
1497	!
1498	!-- Interpolation in time of LSF_DATA for ug, vg and w_subs
1499	nt = 1
1500	DO WHILE ( time > time_vert(nt) )
1501	nt = nt + 1
1502	ENDDO
1503	IF ( time /= time_vert(nt) ) THEN
1504	nt = nt - 1
1505	ENDIF
1506
1507	fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
1508
1509	ug = ug_vert(:,nt) + fac * ( ug_vert(:,nt+1) - ug_vert(:,nt) )
1510	vg = vg_vert(:,nt) + fac * ( vg_vert(:,nt+1) - vg_vert(:,nt) )
1511
1512	IF ( large_scale_subsidence ) THEN
1513	w_subs = wsubs_vert(:,nt) &
1514	+ fac * ( wsubs_vert(:,nt+1) - wsubs_vert(:,nt) )
1515	ENDIF
1516
1517	END SUBROUTINE ls_forcing_vert
1518
1519
1520	!------------------------------------------------------------------------------!
1521	! Description:
1522	! ------------
1523	!> Call for all grid points
1524	!------------------------------------------------------------------------------!
1525	SUBROUTINE ls_advec ( time, prog_var )
1526
1527
1528	IMPLICIT NONE
1529
1530	CHARACTER (LEN=*) :: prog_var !<
1531
1532	REAL(wp), INTENT(in) :: time !<
1533	REAL(wp) :: fac !<
1534
1535	INTEGER(iwp) :: i !<
1536	INTEGER(iwp) :: j !<
1537	INTEGER(iwp) :: k !<
1538	INTEGER(iwp) :: nt !<
1539
1540	!
1541	!-- Interpolation in time of LSF_DATA
1542	nt = 1
1543	DO WHILE ( time > time_vert(nt) )
1544	nt = nt + 1
1545	ENDDO
1546	IF ( time /= time_vert(nt) ) THEN
1547	nt = nt - 1
1548	ENDIF
1549
1550	fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
1551
1552	!
1553	!-- Add horizontal large scale advection tendencies of pt and q
1554	SELECT CASE ( prog_var )
1555
1556	CASE ( 'pt' )
1557
1558	DO i = nxl, nxr
1559	DO j = nys, nyn
1560	DO k = nzb+1, nzt
1561	tend(k,j,i) = tend(k,j,i) + td_lsa_lpt(k,nt) + fac * &
1562	( td_lsa_lpt(k,nt+1) - td_lsa_lpt(k,nt) ) *&
1563	MERGE( 1.0_wp, 0.0_wp, &
1564	BTEST( wall_flags_0(k,j,i), 0 ) )
1565	ENDDO
1566	ENDDO
1567	ENDDO
1568
1569	CASE ( 'q' )
1570
1571	DO i = nxl, nxr
1572	DO j = nys, nyn
1573	DO k = nzb+1, nzt
1574	tend(k,j,i) = tend(k,j,i) + td_lsa_q(k,nt) + fac * &
1575	( td_lsa_q(k,nt+1) - td_lsa_q(k,nt) ) * &
1576	MERGE( 1.0_wp, 0.0_wp, &
1577	BTEST( wall_flags_0(k,j,i), 0 ) )
1578	ENDDO
1579	ENDDO
1580	ENDDO
1581
1582	END SELECT
1583
1584	!
1585	!-- Subsidence of pt and q with prescribed subsidence tendencies
1586	IF ( large_scale_subsidence .AND. use_subsidence_tendencies ) THEN
1587
1588	SELECT CASE ( prog_var )
1589
1590	CASE ( 'pt' )
1591
1592	DO i = nxl, nxr
1593	DO j = nys, nyn
1594	DO k = nzb+1, nzt
1595	tend(k,j,i) = tend(k,j,i) + td_sub_lpt(k,nt) + fac * &
1596	( td_sub_lpt(k,nt+1) - td_sub_lpt(k,nt) )*&
1597	MERGE( 1.0_wp, 0.0_wp, &
1598	BTEST( wall_flags_0(k,j,i), 0 ) )
1599	ENDDO
1600	ENDDO
1601	ENDDO
1602
1603	CASE ( 'q' )
1604
1605	DO i = nxl, nxr
1606	DO j = nys, nyn
1607	DO k = nzb+1, nzt
1608	tend(k,j,i) = tend(k,j,i) + td_sub_q(k,nt) + fac * &
1609	( td_sub_q(k,nt+1) - td_sub_q(k,nt) ) * &
1610	MERGE( 1.0_wp, 0.0_wp, &
1611	BTEST( wall_flags_0(k,j,i), 0 ) )
1612	ENDDO
1613	ENDDO
1614	ENDDO
1615
1616	END SELECT
1617
1618	ENDIF
1619
1620	END SUBROUTINE ls_advec
1621
1622
1623	!------------------------------------------------------------------------------!
1624	! Description:
1625	! ------------
1626	!> Call for grid point i,j
1627	!------------------------------------------------------------------------------!
1628	SUBROUTINE ls_advec_ij ( i, j, time, prog_var )
1629
1630	IMPLICIT NONE
1631
1632	CHARACTER (LEN=*) :: prog_var !<
1633
1634	REAL(wp), INTENT(in) :: time !<
1635	REAL(wp) :: fac !<
1636
1637	INTEGER(iwp) :: i !<
1638	INTEGER(iwp) :: j !<
1639	INTEGER(iwp) :: k !<
1640	INTEGER(iwp) :: nt !<
1641
1642	!
1643	!-- Interpolation in time of LSF_DATA
1644	nt = 1
1645	DO WHILE ( time > time_vert(nt) )
1646	nt = nt + 1
1647	ENDDO
1648	IF ( time /= time_vert(nt) ) THEN
1649	nt = nt - 1
1650	ENDIF
1651
1652	fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
1653
1654	!
1655	!-- Add horizontal large scale advection tendencies of pt and q
1656	SELECT CASE ( prog_var )
1657
1658	CASE ( 'pt' )
1659
1660	DO k = nzb+1, nzt
1661	tend(k,j,i) = tend(k,j,i) + td_lsa_lpt(k,nt) &
1662	+ fac * ( td_lsa_lpt(k,nt+1) - td_lsa_lpt(k,nt) )*&
1663	MERGE( 1.0_wp, 0.0_wp, &
1664	BTEST( wall_flags_0(k,j,i), 0 ) )
1665	ENDDO
1666
1667	CASE ( 'q' )
1668
1669	DO k = nzb+1, nzt
1670	tend(k,j,i) = tend(k,j,i) + td_lsa_q(k,nt) &
1671	+ fac * ( td_lsa_q(k,nt+1) - td_lsa_q(k,nt) ) * &
1672	MERGE( 1.0_wp, 0.0_wp, &
1673	BTEST( wall_flags_0(k,j,i), 0 ) )
1674	ENDDO
1675
1676	END SELECT
1677
1678	!
1679	!-- Subsidence of pt and q with prescribed profiles
1680	IF ( large_scale_subsidence .AND. use_subsidence_tendencies ) THEN
1681
1682	SELECT CASE ( prog_var )
1683
1684	CASE ( 'pt' )
1685
1686	DO k = nzb+1, nzt
1687	tend(k,j,i) = tend(k,j,i) + td_sub_lpt(k,nt) + fac * &
1688	( td_sub_lpt(k,nt+1) - td_sub_lpt(k,nt) ) * &
1689	MERGE( 1.0_wp, 0.0_wp, &
1690	BTEST( wall_flags_0(k,j,i), 0 ) )
1691	ENDDO
1692
1693	CASE ( 'q' )
1694
1695	DO k = nzb+1, nzt
1696	tend(k,j,i) = tend(k,j,i) + td_sub_q(k,nt) + fac * &
1697	( td_sub_q(k,nt+1) - td_sub_q(k,nt) ) * &
1698	MERGE( 1.0_wp, 0.0_wp, &
1699	BTEST( wall_flags_0(k,j,i), 0 ) )
1700	ENDDO
1701
1702	END SELECT
1703
1704	ENDIF
1705
1706	END SUBROUTINE ls_advec_ij
1707
1708
1709	!------------------------------------------------------------------------------!
1710	! Description:
1711	! ------------
1712	!> @todo Missing subroutine description.
1713	!------------------------------------------------------------------------------!
1714	SUBROUTINE nudge_init
1715
1716	IMPLICIT NONE
1717
1718
1719	INTEGER(iwp) :: finput = 90 !<
1720	INTEGER(iwp) :: ierrn !<
1721	INTEGER(iwp) :: k !<
1722	INTEGER(iwp) :: nt !<
1723
1724	CHARACTER(1) :: hash !<
1725
1726	REAL(wp) :: highheight !<
1727	REAL(wp) :: highqnudge !<
1728	REAL(wp) :: highptnudge !<
1729	REAL(wp) :: highunudge !<
1730	REAL(wp) :: highvnudge !<
1731	REAL(wp) :: highwnudge !<
1732	REAL(wp) :: hightnudge !<
1733
1734	REAL(wp) :: lowheight !<
1735	REAL(wp) :: lowqnudge !<
1736	REAL(wp) :: lowptnudge !<
1737	REAL(wp) :: lowunudge !<
1738	REAL(wp) :: lowvnudge !<
1739	REAL(wp) :: lowwnudge !<
1740	REAL(wp) :: lowtnudge !<
1741
1742	REAL(wp) :: fac !<
1743
1744	ALLOCATE( ptnudge(nzb:nzt+1,1:ntnudge), qnudge(nzb:nzt+1,1:ntnudge), &
1745	tnudge(nzb:nzt+1,1:ntnudge), unudge(nzb:nzt+1,1:ntnudge), &
1746	vnudge(nzb:nzt+1,1:ntnudge), wnudge(nzb:nzt+1,1:ntnudge) )
1747
1748	ALLOCATE( tmp_tnudge(nzb:nzt) )
1749
1750	ALLOCATE( timenudge(0:ntnudge) )
1751
1752	ptnudge = 0.0_wp; qnudge = 0.0_wp; tnudge = 0.0_wp; unudge = 0.0_wp
1753	vnudge = 0.0_wp; wnudge = 0.0_wp; timenudge = 0.0_wp
1754	!
1755	!-- Initialize array tmp_nudge with a current nudging time scale of 6 hours
1756	tmp_tnudge = 21600.0_wp
1757
1758	nt = 0
1759	OPEN ( finput, FILE='NUDGING_DATA', STATUS='OLD', &
1760	FORM='FORMATTED', IOSTAT=ierrn )
1761
1762	IF ( ierrn /= 0 ) THEN
1763	message_string = 'file NUDGING_DATA does not exist'
1764	CALL message( 'nudging', 'PA0365', 1, 2, 0, 6, 0 )
1765	ENDIF
1766
1767	ierrn = 0
1768
1769	rloop:DO
1770	nt = nt + 1
1771	hash = "#"
1772	ierrn = 1 ! not zero
1773	!
1774	!-- Search for the next line consisting of "# time",
1775	!-- from there onwards the profiles will be read
1776	DO WHILE ( .NOT. ( hash == "#" .AND. ierrn == 0 ) )
1777
1778	READ ( finput, *, IOSTAT=ierrn ) hash, timenudge(nt)
1779	IF ( ierrn < 0 ) EXIT rloop
1780
1781	ENDDO
1782
1783	ierrn = 0
1784	READ ( finput, *, IOSTAT=ierrn ) lowheight, lowtnudge, lowunudge, &
1785	lowvnudge, lowwnudge , lowptnudge, &
1786	lowqnudge
1787
1788	IF ( ierrn /= 0 ) THEN
1789	message_string = 'errors in file NUDGING_DATA'
1790	CALL message( 'nudging', 'PA0366', 1, 2, 0, 6, 0 )
1791	ENDIF
1792
1793	ierrn = 0
1794	READ ( finput, *, IOSTAT=ierrn ) highheight, hightnudge, highunudge, &
1795	highvnudge, highwnudge , highptnudge, &
1796	highqnudge
1797
1798	IF ( ierrn /= 0 ) THEN
1799	message_string = 'errors in file NUDGING_DATA'
1800	CALL message( 'nudging', 'PA0366', 1, 2, 0, 6, 0 )
1801	ENDIF
1802
1803	DO k = nzb, nzt+1
1804	DO WHILE ( highheight < zu(k) )
1805	lowheight = highheight
1806	lowtnudge = hightnudge
1807	lowunudge = highunudge
1808	lowvnudge = highvnudge
1809	lowwnudge = highwnudge
1810	lowptnudge = highptnudge
1811	lowqnudge = highqnudge
1812
1813	ierrn = 0
1814	READ ( finput, *, IOSTAT=ierrn ) highheight , hightnudge , &
1815	highunudge , highvnudge , &
1816	highwnudge , highptnudge, &
1817	highqnudge
1818	IF (ierrn /= 0 ) THEN
1819	WRITE( message_string, * ) 'zu(',k,') = ', zu(k), 'm is ',&
1820	'higher than the maximum height in NUDING_DATA which ', &
1821	'is ', lowheight, 'm. Interpolation on PALM ', &
1822	'grid is not possible.'
1823	CALL message( 'nudging', 'PA0364', 1, 2, 0, 6, 0 )
1824	ENDIF
1825	ENDDO
1826
1827	!
1828	!-- Interpolation of prescribed profiles in space
1829
1830	fac = ( highheight - zu(k) ) / ( highheight - lowheight )
1831
1832	tnudge(k,nt) = fac * lowtnudge + ( 1.0_wp - fac ) * hightnudge
1833	unudge(k,nt) = fac * lowunudge + ( 1.0_wp - fac ) * highunudge
1834	vnudge(k,nt) = fac * lowvnudge + ( 1.0_wp - fac ) * highvnudge
1835	wnudge(k,nt) = fac * lowwnudge + ( 1.0_wp - fac ) * highwnudge
1836	ptnudge(k,nt) = fac * lowptnudge + ( 1.0_wp - fac ) * highptnudge
1837	qnudge(k,nt) = fac * lowqnudge + ( 1.0_wp - fac ) * highqnudge
1838	ENDDO
1839
1840	ENDDO rloop
1841
1842	CLOSE ( finput )
1843
1844	!
1845	!-- Overwrite initial profiles in case of nudging
1846	IF ( nudging ) THEN
1847	pt_init = ptnudge(:,1)
1848	u_init = unudge(:,1)
1849	v_init = vnudge(:,1)
1850	IF ( humidity ) THEN ! is passive_scalar correct???
1851	q_init = qnudge(:,1)
1852	ENDIF
1853
1854	WRITE( message_string, * ) 'Initial profiles of u, v and ', &
1855	'scalars from NUDGING_DATA are used.'
1856	CALL message( 'large_scale_forcing_nudging', 'PA0370', 0, 0, 0, 6, 0 )
1857	ENDIF
1858
1859
1860	END SUBROUTINE nudge_init
1861
1862	!------------------------------------------------------------------------------!
1863	! Description:
1864	! ------------
1865	!> @todo Missing subroutine description.
1866	!------------------------------------------------------------------------------!
1867	SUBROUTINE calc_tnudge ( time )
1868
1869	IMPLICIT NONE
1870
1871
1872	REAL(wp) :: dtm !<
1873	REAL(wp) :: dtp !<
1874	REAL(wp) :: time !<
1875
1876	INTEGER(iwp) :: k !<
1877	INTEGER(iwp) :: nt !<
1878
1879	nt = 1
1880	DO WHILE ( time > timenudge(nt) )
1881	nt = nt+1
1882	ENDDO
1883	IF ( time /= timenudge(1) ) THEN
1884	nt = nt-1
1885	ENDIF
1886
1887	dtm = ( time - timenudge(nt) ) / ( timenudge(nt+1) - timenudge(nt) )
1888	dtp = ( timenudge(nt+1) - time ) / ( timenudge(nt+1) - timenudge(nt) )
1889
1890	DO k = nzb, nzt
1891	tmp_tnudge(k) = MAX( dt_3d, tnudge(k,nt) * dtp + tnudge(k,nt+1) * dtm )
1892	ENDDO
1893
1894	END SUBROUTINE calc_tnudge
1895
1896	!------------------------------------------------------------------------------!
1897	! Description:
1898	! ------------
1899	!> Call for all grid points
1900	!------------------------------------------------------------------------------!
1901	SUBROUTINE nudge ( time, prog_var )
1902
1903	IMPLICIT NONE
1904
1905	CHARACTER (LEN=*) :: prog_var !<
1906
1907	REAL(wp) :: tmp_tend !<
1908	REAL(wp) :: dtm !<
1909	REAL(wp) :: dtp !<
1910	REAL(wp) :: time !<
1911
1912	INTEGER(iwp) :: i !<
1913	INTEGER(iwp) :: j !<
1914	INTEGER(iwp) :: k !<
1915	INTEGER(iwp) :: nt !<
1916
1917
1918	nt = 1
1919	DO WHILE ( time > timenudge(nt) )
1920	nt = nt+1
1921	ENDDO
1922	IF ( time /= timenudge(1) ) THEN
1923	nt = nt-1
1924	ENDIF
1925
1926	dtm = ( time - timenudge(nt) ) / ( timenudge(nt+1) - timenudge(nt) )
1927	dtp = ( timenudge(nt+1) - time ) / ( timenudge(nt+1) - timenudge(nt) )
1928
1929	SELECT CASE ( prog_var )
1930
1931	CASE ( 'u' )
1932
1933	DO i = nxl, nxr
1934	DO j = nys, nyn
1935
1936	DO k = nzb+1, nzt
1937
1938	tmp_tend = - ( hom(k,1,1,0) - ( unudge(k,nt) * dtp + &
1939	unudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
1940
1941	tend(k,j,i) = tend(k,j,i) + tmp_tend * &
1942	MERGE( 1.0_wp, 0.0_wp, &
1943	BTEST( wall_flags_0(k,j,i), 1 ) )
1944
1945	sums_ls_l(k,6) = sums_ls_l(k,6) + tmp_tend * &
1946	weight_substep(intermediate_timestep_count)
1947	ENDDO
1948
1949	sums_ls_l(nzt+1,6) = sums_ls_l(nzt,6)
1950
1951	ENDDO
1952	ENDDO
1953
1954	CASE ( 'v' )
1955
1956	DO i = nxl, nxr
1957	DO j = nys, nyn
1958
1959	DO k = nzb+1, nzt
1960
1961	tmp_tend = - ( hom(k,1,2,0) - ( vnudge(k,nt) * dtp + &
1962	vnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
1963
1964	tend(k,j,i) = tend(k,j,i) + tmp_tend * &
1965	MERGE( 1.0_wp, 0.0_wp, &
1966	BTEST( wall_flags_0(k,j,i), 2 ) )
1967
1968	sums_ls_l(k,7) = sums_ls_l(k,7) + tmp_tend * &
1969	weight_substep(intermediate_timestep_count)
1970	ENDDO
1971
1972	sums_ls_l(nzt+1,7) = sums_ls_l(nzt,7)
1973
1974	ENDDO
1975	ENDDO
1976
1977	CASE ( 'pt' )
1978
1979	DO i = nxl, nxr
1980	DO j = nys, nyn
1981
1982	DO k = nzb+1, nzt
1983
1984	tmp_tend = - ( hom(k,1,4,0) - ( ptnudge(k,nt) * dtp + &
1985	ptnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
1986
1987	tend(k,j,i) = tend(k,j,i) + tmp_tend * &
1988	MERGE( 1.0_wp, 0.0_wp, &
1989	BTEST( wall_flags_0(k,j,i), 0 ) )
1990
1991	sums_ls_l(k,4) = sums_ls_l(k,4) + tmp_tend * &
1992	weight_substep(intermediate_timestep_count)
1993	ENDDO
1994
1995	sums_ls_l(nzt+1,4) = sums_ls_l(nzt,4)
1996
1997	ENDDO
1998	ENDDO
1999
2000	CASE ( 'q' )
2001
2002	DO i = nxl, nxr
2003	DO j = nys, nyn
2004
2005	DO k = nzb+1, nzt
2006
2007	tmp_tend = - ( hom(k,1,41,0) - ( qnudge(k,nt) * dtp + &
2008	qnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
2009
2010	tend(k,j,i) = tend(k,j,i) + tmp_tend * &
2011	MERGE( 1.0_wp, 0.0_wp, &
2012	BTEST( wall_flags_0(k,j,i), 0 ) )
2013
2014	sums_ls_l(k,5) = sums_ls_l(k,5) + tmp_tend * &
2015	weight_substep(intermediate_timestep_count)
2016	ENDDO
2017
2018	sums_ls_l(nzt+1,5) = sums_ls_l(nzt,5)
2019
2020	ENDDO
2021	ENDDO
2022
2023	CASE DEFAULT
2024	message_string = 'unknown prognostic variable "' // prog_var // '"'
2025	CALL message( 'nudge', 'PA0367', 1, 2, 0, 6, 0 )
2026
2027	END SELECT
2028
2029	END SUBROUTINE nudge
2030
2031
2032	!------------------------------------------------------------------------------!
2033	! Description:
2034	! ------------
2035	!> Call for grid point i,j
2036	!------------------------------------------------------------------------------!
2037
2038	SUBROUTINE nudge_ij( i, j, time, prog_var )
2039
2040	IMPLICIT NONE
2041
2042
2043	CHARACTER (LEN=*) :: prog_var !<
2044
2045	REAL(wp) :: tmp_tend !<
2046	REAL(wp) :: dtm !<
2047	REAL(wp) :: dtp !<
2048	REAL(wp) :: time !<
2049
2050	INTEGER(iwp) :: i !<
2051	INTEGER(iwp) :: j !<
2052	INTEGER(iwp) :: k !<
2053	INTEGER(iwp) :: nt !<
2054
2055
2056	nt = 1
2057	DO WHILE ( time > timenudge(nt) )
2058	nt = nt+1
2059	ENDDO
2060	IF ( time /= timenudge(1) ) THEN
2061	nt = nt-1
2062	ENDIF
2063
2064	dtm = ( time - timenudge(nt) ) / ( timenudge(nt+1) - timenudge(nt) )
2065	dtp = ( timenudge(nt+1) - time ) / ( timenudge(nt+1) - timenudge(nt) )
2066
2067	SELECT CASE ( prog_var )
2068
2069	CASE ( 'u' )
2070
2071	DO k = nzb+1, nzt
2072
2073	tmp_tend = - ( hom(k,1,1,0) - ( unudge(k,nt) * dtp + &
2074	unudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
2075
2076	tend(k,j,i) = tend(k,j,i) + tmp_tend * &
2077	MERGE( 1.0_wp, 0.0_wp, &
2078	BTEST( wall_flags_0(k,j,i), 1 ) )
2079
2080	sums_ls_l(k,6) = sums_ls_l(k,6) + tmp_tend &
2081	* weight_substep(intermediate_timestep_count)
2082	ENDDO
2083
2084	sums_ls_l(nzt+1,6) = sums_ls_l(nzt,6)
2085
2086	CASE ( 'v' )
2087
2088	DO k = nzb+1, nzt
2089
2090	tmp_tend = - ( hom(k,1,2,0) - ( vnudge(k,nt) * dtp + &
2091	vnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
2092
2093	tend(k,j,i) = tend(k,j,i) + tmp_tend * &
2094	MERGE( 1.0_wp, 0.0_wp, &
2095	BTEST( wall_flags_0(k,j,i), 2 ) )
2096
2097	sums_ls_l(k,7) = sums_ls_l(k,7) + tmp_tend &
2098	* weight_substep(intermediate_timestep_count)
2099	ENDDO
2100
2101	sums_ls_l(nzt+1,7) = sums_ls_l(nzt,7)
2102
2103	CASE ( 'pt' )
2104
2105	DO k = nzb+1, nzt
2106
2107	tmp_tend = - ( hom(k,1,4,0) - ( ptnudge(k,nt) * dtp + &
2108	ptnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
2109
2110	tend(k,j,i) = tend(k,j,i) + tmp_tend * &
2111	MERGE( 1.0_wp, 0.0_wp, &
2112	BTEST( wall_flags_0(k,j,i), 0 ) )
2113
2114	sums_ls_l(k,4) = sums_ls_l(k,4) + tmp_tend &
2115	* weight_substep(intermediate_timestep_count)
2116	ENDDO
2117
2118	sums_ls_l(nzt+1,4) = sums_ls_l(nzt,4)
2119
2120
2121	CASE ( 'q' )
2122
2123	DO k = nzb+1, nzt
2124
2125	tmp_tend = - ( hom(k,1,41,0) - ( qnudge(k,nt) * dtp + &
2126	qnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
2127
2128	tend(k,j,i) = tend(k,j,i) + tmp_tend * &
2129	MERGE( 1.0_wp, 0.0_wp, &
2130	BTEST( wall_flags_0(k,j,i), 0 ) )
2131
2132	sums_ls_l(k,5) = sums_ls_l(k,5) + tmp_tend &
2133	* weight_substep(intermediate_timestep_count)
2134	ENDDO
2135
2136	sums_ls_l(nzt+1,5) = sums_ls_l(nzt,5)
2137
2138	CASE DEFAULT
2139	message_string = 'unknown prognostic variable "' // prog_var // '"'
2140	CALL message( 'nudge', 'PA0367', 1, 2, 0, 6, 0 )
2141
2142	END SELECT
2143
2144
2145	END SUBROUTINE nudge_ij
2146
2147
2148	!------------------------------------------------------------------------------!
2149	! Description:
2150	! ------------
2151	!> @todo Missing subroutine description.
2152	!------------------------------------------------------------------------------!
2153	SUBROUTINE nudge_ref ( time )
2154
2155	IMPLICIT NONE
2156
2157	INTEGER(iwp) :: nt !<
2158
2159	REAL(wp) :: fac !<
2160	REAL(wp), INTENT(in) :: time !<
2161
2162	!
2163	!-- Interpolation in time of NUDGING_DATA for pt_init and q_init. This is
2164	!-- needed for correct upper boundary conditions for pt and q and in case that
2165	! large scale subsidence as well as scalar Rayleigh-damping are used
2166	nt = 1
2167	DO WHILE ( time > time_vert(nt) )
2168	nt = nt + 1
2169	ENDDO
2170	IF ( time /= time_vert(nt) ) THEN
2171	nt = nt - 1
2172	ENDIF
2173
2174	fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
2175
2176	pt_init = ptnudge(:,nt) + fac * ( ptnudge(:,nt+1) - ptnudge(:,nt) )
2177	q_init = qnudge(:,nt) + fac * ( qnudge(:,nt+1) - qnudge(:,nt) )
2178	u_init = unudge(:,nt) + fac * ( unudge(:,nt+1) - unudge(:,nt) )
2179	v_init = vnudge(:,nt) + fac * ( vnudge(:,nt+1) - vnudge(:,nt) )
2180
2181	END SUBROUTINE nudge_ref
2182
2183
2184	END MODULE lsf_nudging_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |