Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |