Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |