Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |