Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |