Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |