Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |