Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |