Home

Context Navigation

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

Last change on this file since 2347 was 2342, checked in by boeske, 7 years ago
fixed check if surface forcing data is available until end of simulation
Property svn:keywords set to `Id`
File size: 53.0 KB

Line
1	!> @file lsf_nudging_mod.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of PALM.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the
6	! terms of the GNU General Public License as published by the Free Software
7	! Foundation, either version 3 of the License, or (at your option) any later
8	! version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 1997-2017 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! ------------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: large_scale_forcing_nudging_mod.f90 2342 2017-08-08 11:00:43Z knoop $
27	! fixed check if surface forcing data is available until end of simulation
28	!
29	! 2320 2017-07-21 12:47:43Z suehring
30	! initial revision
31	!
32	! Description:
33	! ------------
34	!> Calculates large scale forcings (geostrophic wind and subsidence velocity) as
35	!> well as surfaces fluxes dependent on time given in an external file (LSF_DATA).
36	!> Moreover, module contains nudging routines, where u, v, pt and q are nudged
37	!> to given profiles on a relaxation timescale tnudge.
38	!> Profiles are read in from NUDGING_DATA.
39	!> Code is based on Neggers et al. (2012) and also in parts on DALES and UCLA-LES.
40	!> @todo: Revise reading of ASCII-files
41	!> @todo: Remove unused variables and control flags
42	!> @todo: Revise large-scale facing of surface variables
43	!> @todo: Revise control flags lsf_exception, lsf_surf, lsf_vert, etc.
44	!--------------------------------------------------------------------------------!
45	MODULE lsf_nudging_mod
46
47	USE arrays_3d, &
48	ONLY: heatflux_input_conversion, pt, pt_init, q, q_init, tend, &
49	u, u_init, ug, v, v_init, vg, w_subs, &
50	waterflux_input_conversion, zu
51
52	USE control_parameters, &
53	ONLY: bc_lr, bc_ns, bc_pt_b, bc_q_b, data_output_pr, dt_3d, end_time, &
54	humidity, intermediate_timestep_count, ibc_pt_b, ibc_q_b, &
55	large_scale_forcing, large_scale_subsidence, lsf_surf, lsf_vert,&
56	lsf_exception, message_string, nudging, passive_scalar, &
57	pt_surface, ocean, q_surface, surface_pressure, topography, &
58	use_subsidence_tendencies
59
60
61	USE indices, &
62	ONLY: ngp_sums_ls, nxl, nxr, nys, nyn, nzb, nz, nzt, wall_flags_0
63
64	USE kinds
65
66	USE surface_mod, &
67	ONLY: surf_def_h, surf_lsm_h, surf_usm_h
68
69	USE statistics, &
70	ONLY: hom, statistic_regions, sums_ls_l, weight_substep
71
72
73	INTEGER(iwp) :: nlsf = 1000 !< maximum number of profiles in LSF_DATA (large scale forcing)
74	INTEGER(iwp) :: ntnudge = 1000 !< maximum number of profiles in NUDGING_DATA (nudging)
75
76	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ptnudge !< vertical profile of pot. temperature interpolated to vertical grid (nudging)
77	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: qnudge !< vertical profile of specific humidity interpolated to vertical grid (nudging)
78	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tnudge !< vertical profile of nudging time scale interpolated to vertical grid (nudging)
79	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: td_lsa_lpt !< temperature tendency due to large scale advection (large scale forcing)
80	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: td_lsa_q !< specific humidity tendency due to large scale advection (large scale forcing)
81	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: td_sub_lpt !< temperature tendency due to subsidence/ascent (large scale forcing)
82	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: td_sub_q !< specific humidity tendency due to subsidence/ascent (large scale forcing)
83	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ug_vert !< vertical profile of geostrophic wind component in x-direction interpolated to vertical grid (large scale forcing)
84	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: unudge !< vertical profile of wind component in x-direction interpolated to vertical grid (nudging)
85	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: vnudge !< vertical profile of wind component in y-direction interpolated to vertical grid (nudging)
86	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: vg_vert !< vertical profile of geostrophic wind component in y-direction interpolated to vertical grid (large scale forcing)
87	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: wnudge !< vertical profile of subsidence/ascent velocity interpolated to vertical grid (nudging) ???
88	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: wsubs_vert !< vertical profile of wind component in z-direction interpolated to vertical grid (nudging) ???
89
90	REAL(wp), DIMENSION(:), ALLOCATABLE :: shf_surf !< time-dependent surface sensible heat flux (large scale forcing)
91	REAL(wp), DIMENSION(:), ALLOCATABLE :: timenudge !< times at which vertical profiles are defined in NUDGING_DATA (nudging)
92	REAL(wp), DIMENSION(:), ALLOCATABLE :: time_surf !< times at which surface values/fluxes are defined in LSF_DATA (large scale forcing)
93	REAL(wp), DIMENSION(:), ALLOCATABLE :: time_vert !< times at which vertical profiles are defined in LSF_DATA (large scale forcing)
94
95	REAL(wp), DIMENSION(:), ALLOCATABLE :: tmp_tnudge !< current nudging time scale
96
97	REAL(wp), DIMENSION(:), ALLOCATABLE :: p_surf !< time-dependent surface pressure (large scale forcing)
98	REAL(wp), DIMENSION(:), ALLOCATABLE :: pt_surf !< time-dependent surface temperature (large scale forcing)
99	REAL(wp), DIMENSION(:), ALLOCATABLE :: qsws_surf !< time-dependent surface latent heat flux (large scale forcing)
100	REAL(wp), DIMENSION(:), ALLOCATABLE :: q_surf !< time-dependent surface specific humidity (large scale forcing)
101
102	SAVE
103	PRIVATE
104	!
105	!-- Public subroutines
106	PUBLIC ls_forcing_surf, ls_forcing_vert, ls_advec, lsf_init, &
107	lsf_nudging_check_parameters, nudge_init, &
108	lsf_nudging_check_data_output_pr, lsf_nudging_header, &
109	calc_tnudge, nudge, nudge_ref
110	!
111	!-- Public variables
112	PUBLIC qsws_surf, shf_surf, td_lsa_lpt, td_lsa_q, td_sub_lpt, &
113	td_sub_q, time_vert
114
115	INTERFACE ls_advec
116	MODULE PROCEDURE ls_advec
117	MODULE PROCEDURE ls_advec_ij
118	END INTERFACE ls_advec
119
120	INTERFACE nudge
121	MODULE PROCEDURE nudge
122	MODULE PROCEDURE nudge_ij
123	END INTERFACE nudge
124
125	CONTAINS
126
127
128	!------------------------------------------------------------------------------!
129	! Description:
130	! ------------
131	!> @todo Missing subroutine description.
132	!------------------------------------------------------------------------------!
133	SUBROUTINE lsf_nudging_check_parameters
134
135	IMPLICIT NONE
136	!
137	!-- Check nudging and large scale forcing from external file
138	IF ( nudging .AND. ( .NOT. large_scale_forcing ) ) THEN
139	message_string = 'Nudging requires large_scale_forcing = .T.. &'// &
140	'Surface fluxes and geostrophic wind should be &'// &
141	'prescribed in file LSF_DATA'
142	CALL message( 'check_parameters', 'PA0374', 1, 2, 0, 6, 0 )
143	ENDIF
144
145	IF ( large_scale_forcing .AND. ( bc_lr /= 'cyclic' .OR. &
146	bc_ns /= 'cyclic' ) ) THEN
147	message_string = 'Non-cyclic lateral boundaries do not allow for &'//&
148	'the usage of large scale forcing from external file.'
149	CALL message( 'check_parameters', 'PA0375', 1, 2, 0, 6, 0 )
150	ENDIF
151
152	IF ( large_scale_forcing .AND. ( .NOT. humidity ) ) THEN
153	message_string = 'The usage of large scale forcing from external &'//&
154	'file LSF_DATA requires humidity = .T..'
155	CALL message( 'check_parameters', 'PA0376', 1, 2, 0, 6, 0 )
156	ENDIF
157
158	IF ( large_scale_forcing .AND. passive_scalar ) THEN
159	message_string = 'The usage of large scale forcing from external &'//&
160	'file LSF_DATA is not implemented for passive scalars'
161	CALL message( 'check_parameters', 'PA0440', 1, 2, 0, 6, 0 )
162	ENDIF
163
164	IF ( large_scale_forcing .AND. topography /= 'flat' &
165	.AND. .NOT. lsf_exception ) THEN
166	message_string = 'The usage of large scale forcing from external &'//&
167	'file LSF_DATA is not implemented for non-flat topography'
168	CALL message( 'check_parameters', 'PA0377', 1, 2, 0, 6, 0 )
169	ENDIF
170
171	IF ( large_scale_forcing .AND. ocean ) THEN
172	message_string = 'The usage of large scale forcing from external &'//&
173	'file LSF_DATA is not implemented for ocean runs'
174	CALL message( 'check_parameters', 'PA0378', 1, 2, 0, 6, 0 )
175	ENDIF
176
177	END SUBROUTINE lsf_nudging_check_parameters
178
179	!------------------------------------------------------------------------------!
180	! Description:
181	! ------------
182	!> Check data output of profiles for land surface model
183	!------------------------------------------------------------------------------!
184	SUBROUTINE lsf_nudging_check_data_output_pr( variable, var_count, unit, &
185	dopr_unit )
186
187	USE profil_parameter
188
189	IMPLICIT NONE
190
191	CHARACTER (LEN=*) :: unit !<
192	CHARACTER (LEN=*) :: variable !<
193	CHARACTER (LEN=*) :: dopr_unit !< local value of dopr_unit
194
195	INTEGER(iwp) :: user_pr_index !<
196	INTEGER(iwp) :: var_count !<
197
198	SELECT CASE ( TRIM( variable ) )
199
200
201	CASE ( 'td_lsa_lpt' )
202	IF ( .NOT. large_scale_forcing ) THEN
203	message_string = 'data_output_pr = ' // &
204	TRIM( data_output_pr(var_count) ) // &
205	' is not implemented for ' // &
206	'large_scale_forcing = .FALSE.'
207	CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393', &
208	1, 2, 0, 6, 0 )
209	ELSE
210	dopr_index(var_count) = 81
211	dopr_unit = 'K/s'
212	unit = 'K/s'
213	hom(:,2,81,:) = SPREAD( zu, 2, statistic_regions+1 )
214	ENDIF
215
216	CASE ( 'td_lsa_q' )
217	IF ( .NOT. large_scale_forcing ) THEN
218	message_string = 'data_output_pr = ' // &
219	TRIM( data_output_pr(var_count) ) // &
220	' is not implemented for ' // &
221	'large_scale_forcing = .FALSE.'
222	CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393', &
223	1, 2, 0, 6, 0 )
224	ELSE
225	dopr_index(var_count) = 82
226	dopr_unit = 'kg/kgs'
227	unit = 'kg/kgs'
228	hom(:,2,82,:) = SPREAD( zu, 2, statistic_regions+1 )
229	ENDIF
230	CASE ( 'td_sub_lpt' )
231	IF ( .NOT. large_scale_forcing ) THEN
232	message_string = 'data_output_pr = ' // &
233	TRIM( data_output_pr(var_count) ) // &
234	' is not implemented for ' // &
235	'large_scale_forcing = .FALSE.'
236	CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393', &
237	1, 2, 0, 6, 0 )
238	ELSE
239	dopr_index(var_count) = 83
240	dopr_unit = 'K/s'
241	unit = 'K/s'
242	hom(:,2,83,:) = SPREAD( zu, 2, statistic_regions+1 )
243	ENDIF
244
245	CASE ( 'td_sub_q' )
246	IF ( .NOT. large_scale_forcing ) THEN
247	message_string = 'data_output_pr = ' // &
248	TRIM( data_output_pr(var_count) ) // &
249	' is not implemented for ' // &
250	'large_scale_forcing = .FALSE.'
251	CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393', &
252	1, 2, 0, 6, 0 )
253	ELSE
254	dopr_index(var_count) = 84
255	dopr_unit = 'kg/kgs'
256	unit = 'kg/kgs'
257	hom(:,2,84,:) = SPREAD( zu, 2, statistic_regions+1 )
258	ENDIF
259
260	CASE ( 'td_nud_lpt' )
261	IF ( .NOT. nudging ) THEN
262	message_string = 'data_output_pr = ' // &
263	TRIM( data_output_pr(var_count) ) // &
264	' is not implemented for ' // &
265	'large_scale_forcing = .FALSE.'
266	CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393', &
267	1, 2, 0, 6, 0 )
268	ELSE
269	dopr_index(var_count) = 85
270	dopr_unit = 'K/s'
271	unit = 'K/s'
272	hom(:,2,85,:) = SPREAD( zu, 2, statistic_regions+1 )
273	ENDIF
274
275	CASE ( 'td_nud_q' )
276	IF ( .NOT. nudging ) THEN
277	message_string = 'data_output_pr = ' // &
278	TRIM( data_output_pr(var_count) ) // &
279	' is not implemented for ' // &
280	'large_scale_forcing = .FALSE.'
281	CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393', &
282	1, 2, 0, 6, 0 )
283	ELSE
284	dopr_index(var_count) = 86
285	dopr_unit = 'kg/kgs'
286	unit = 'kg/kgs'
287	hom(:,2,86,:) = SPREAD( zu, 2, statistic_regions+1 )
288	ENDIF
289
290	CASE ( 'td_nud_u' )
291	IF ( .NOT. nudging ) THEN
292	message_string = 'data_output_pr = ' // &
293	TRIM( data_output_pr(var_count) ) // &
294	' is not implemented for ' // &
295	'large_scale_forcing = .FALSE.'
296	CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393', &
297	1, 2, 0, 6, 0 )
298	ELSE
299	dopr_index(var_count) = 87
300	dopr_unit = 'm/s2'
301	unit = 'm/s2'
302	hom(:,2,87,:) = SPREAD( zu, 2, statistic_regions+1 )
303	ENDIF
304
305	CASE ( 'td_nud_v' )
306	IF ( .NOT. nudging ) THEN
307	message_string = 'data_output_pr = ' // &
308	TRIM( data_output_pr(var_count) ) // &
309	' is not implemented for ' // &
310	'large_scale_forcing = .FALSE.'
311	CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393', &
312	1, 2, 0, 6, 0 )
313	ELSE
314	dopr_index(var_count) = 88
315	dopr_unit = 'm/s2'
316	unit = 'm/s2'
317	hom(:,2,88,:) = SPREAD( zu, 2, statistic_regions+1 )
318	ENDIF
319
320
321	CASE DEFAULT
322	unit = 'illegal'
323
324	END SELECT
325
326	END SUBROUTINE lsf_nudging_check_data_output_pr
327
328	!------------------------------------------------------------------------------!
329	! Description:
330	! ------------
331	!> @todo Missing subroutine description.
332	!------------------------------------------------------------------------------!
333	SUBROUTINE lsf_nudging_header ( io )
334
335	IMPLICIT NONE
336
337	INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file
338
339	WRITE ( io, 1 )
340	IF ( large_scale_forcing ) THEN
341	WRITE ( io, 3 )
342	WRITE ( io, 4 )
343
344	IF ( large_scale_subsidence ) THEN
345	IF ( .NOT. use_subsidence_tendencies ) THEN
346	WRITE ( io, 5 )
347	ELSE
348	WRITE ( io, 6 )
349	ENDIF
350	ENDIF
351
352	IF ( bc_pt_b == 'dirichlet' ) THEN
353	WRITE ( io, 12 )
354	ELSEIF ( bc_pt_b == 'neumann' ) THEN
355	WRITE ( io, 13 )
356	ENDIF
357
358	IF ( bc_q_b == 'dirichlet' ) THEN
359	WRITE ( io, 14 )
360	ELSEIF ( bc_q_b == 'neumann' ) THEN
361	WRITE ( io, 15 )
362	ENDIF
363
364	WRITE ( io, 7 )
365	IF ( nudging ) THEN
366	WRITE ( io, 10 )
367	ENDIF
368	ELSE
369	WRITE ( io, 2 )
370	WRITE ( io, 11 )
371	ENDIF
372	IF ( large_scale_subsidence ) THEN
373	WRITE ( io, 8 )
374	WRITE ( io, 9 )
375	ENDIF
376
377
378	1 FORMAT (//' Large scale forcing and nudging:'/ &
379	' -------------------------------'/)
380	2 FORMAT (' --> No large scale forcing from external is used (default) ')
381	3 FORMAT (' --> Large scale forcing from external file LSF_DATA is used: ')
382	4 FORMAT (' - large scale advection tendencies ')
383	5 FORMAT (' - large scale subsidence velocity w_subs ')
384	6 FORMAT (' - large scale subsidence tendencies ')
385	7 FORMAT (' - and geostrophic wind components ug and vg')
386	8 FORMAT (' --> Large-scale vertical motion is used in the ', &
387	'prognostic equation(s) for')
388	9 FORMAT (' the scalar(s) only')
389	10 FORMAT (' --> Nudging is used')
390	11 FORMAT (' --> No nudging is used (default) ')
391	12 FORMAT (' - prescribed surface values for temperature')
392	13 FORMAT (' - prescribed surface fluxes for temperature')
393	14 FORMAT (' - prescribed surface values for humidity')
394	15 FORMAT (' - prescribed surface fluxes for humidity')
395
396	END SUBROUTINE lsf_nudging_header
397
398	!------------------------------------------------------------------------------!
399	! Description:
400	! ------------
401	!> @todo Missing subroutine description.
402	!------------------------------------------------------------------------------!
403	SUBROUTINE lsf_init
404
405	IMPLICIT NONE
406
407	CHARACTER(100) :: chmess !<
408	CHARACTER(1) :: hash !<
409
410	INTEGER(iwp) :: ierrn !<
411	INTEGER(iwp) :: finput = 90 !<
412	INTEGER(iwp) :: k !<
413	INTEGER(iwp) :: nt !<
414
415	REAL(wp) :: fac !<
416	REAL(wp) :: highheight !<
417	REAL(wp) :: highug_vert !<
418	REAL(wp) :: highvg_vert !<
419	REAL(wp) :: highwsubs_vert !<
420	REAL(wp) :: lowheight !<
421	REAL(wp) :: lowug_vert !<
422	REAL(wp) :: lowvg_vert !<
423	REAL(wp) :: lowwsubs_vert !<
424	REAL(wp) :: high_td_lsa_lpt !<
425	REAL(wp) :: low_td_lsa_lpt !<
426	REAL(wp) :: high_td_lsa_q !<
427	REAL(wp) :: low_td_lsa_q !<
428	REAL(wp) :: high_td_sub_lpt !<
429	REAL(wp) :: low_td_sub_lpt !<
430	REAL(wp) :: high_td_sub_q !<
431	REAL(wp) :: low_td_sub_q !<
432	REAL(wp) :: r_dummy !<
433
434	ALLOCATE( p_surf(0:nlsf), pt_surf(0:nlsf), q_surf(0:nlsf), &
435	qsws_surf(0:nlsf), shf_surf(0:nlsf), &
436	td_lsa_lpt(nzb:nzt+1,0:nlsf), td_lsa_q(nzb:nzt+1,0:nlsf), &
437	td_sub_lpt(nzb:nzt+1,0:nlsf), td_sub_q(nzb:nzt+1,0:nlsf), &
438	time_vert(0:nlsf), time_surf(0:nlsf), &
439	ug_vert(nzb:nzt+1,0:nlsf), vg_vert(nzb:nzt+1,0:nlsf), &
440	wsubs_vert(nzb:nzt+1,0:nlsf) )
441
442	p_surf = 0.0_wp; pt_surf = 0.0_wp; q_surf = 0.0_wp; qsws_surf = 0.0_wp
443	shf_surf = 0.0_wp; time_vert = 0.0_wp; td_lsa_lpt = 0.0_wp
444	td_lsa_q = 0.0_wp; td_sub_lpt = 0.0_wp; td_sub_q = 0.0_wp
445	time_surf = 0.0_wp; ug_vert = 0.0_wp; vg_vert = 0.0_wp
446	wsubs_vert = 0.0_wp
447
448	!
449	!-- Array for storing large scale forcing and nudging tendencies at each
450	!-- timestep for data output
451	ALLOCATE( sums_ls_l(nzb:nzt+1,0:7) )
452	sums_ls_l = 0.0_wp
453
454	ngp_sums_ls = (nz+2)*6
455
456	OPEN ( finput, FILE='LSF_DATA', STATUS='OLD', &
457	FORM='FORMATTED', IOSTAT=ierrn )
458
459	IF ( ierrn /= 0 ) THEN
460	message_string = 'file LSF_DATA does not exist'
461	CALL message( 'ls_forcing', 'PA0368', 1, 2, 0, 6, 0 )
462	ENDIF
463
464	ierrn = 0
465	!
466	!-- First three lines of LSF_DATA contain header
467	READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
468	READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
469	READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
470
471	IF ( ierrn /= 0 ) THEN
472	message_string = 'errors in file LSF_DATA'
473	CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
474	ENDIF
475
476	!
477	!-- Surface values are read in
478	nt = 0
479	ierrn = 0
480
481	DO WHILE ( time_surf(nt) < end_time )
482	nt = nt + 1
483	READ ( finput, *, IOSTAT = ierrn ) time_surf(nt), shf_surf(nt), &
484	qsws_surf(nt), pt_surf(nt), &
485	q_surf(nt), p_surf(nt)
486
487	IF ( ierrn /= 0 ) THEN
488	WRITE ( message_string, * ) 'No time dependent surface variables ',&
489	'in&LSF_DATA for end of run found'
490
491	CALL message( 'ls_forcing', 'PA0363', 1, 2, 0, 6, 0 )
492	ENDIF
493	ENDDO
494
495	IF ( time_surf(1) > end_time ) THEN
496	WRITE ( message_string, * ) 'Time dependent surface variables in ', &
497	'&LSF_DATA set in after end of ' , &
498	'simulation - lsf_surf is set to FALSE'
499	CALL message( 'ls_forcing', 'PA0371', 0, 0, 0, 6, 0 )
500	lsf_surf = .FALSE.
501	ENDIF
502
503	!
504	!-- Go to the end of the list with surface variables
505	DO WHILE ( ierrn == 0 )
506	READ ( finput, *, IOSTAT = ierrn ) r_dummy
507	ENDDO
508
509	!
510	!-- Profiles of ug, vg and w_subs are read in (large scale forcing)
511
512	nt = 0
513	DO WHILE ( time_vert(nt) < end_time )
514	nt = nt + 1
515	hash = "#"
516	ierrn = 1 ! not zero
517	!
518	!-- Search for the next line consisting of "# time",
519	!-- from there onwards the profiles will be read
520	DO WHILE ( .NOT. ( hash == "#" .AND. ierrn == 0 ) )
521	READ ( finput, *, IOSTAT=ierrn ) hash, time_vert(nt)
522	IF ( ierrn < 0 ) THEN
523	WRITE( message_string, * ) 'No time dependent vertical profiles',&
524	' in&LSF_DATA for end of run found'
525	CALL message( 'ls_forcing', 'PA0372', 1, 2, 0, 6, 0 )
526	ENDIF
527	ENDDO
528
529	IF ( nt == 1 .AND. time_vert(nt) > end_time ) EXIT
530
531	READ ( finput, *, IOSTAT=ierrn ) lowheight, lowug_vert, lowvg_vert, &
532	lowwsubs_vert, low_td_lsa_lpt, &
533	low_td_lsa_q, low_td_sub_lpt, &
534	low_td_sub_q
535	IF ( ierrn /= 0 ) THEN
536	message_string = 'errors in file LSF_DATA'
537	CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
538	ENDIF
539
540	READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert, &
541	highvg_vert, highwsubs_vert, &
542	high_td_lsa_lpt, high_td_lsa_q, &
543	high_td_sub_lpt, high_td_sub_q
544
545	IF ( ierrn /= 0 ) THEN
546	message_string = 'errors in file LSF_DATA'
547	CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
548	ENDIF
549
550
551	DO k = nzb, nzt+1
552	IF ( highheight < zu(k) ) THEN
553	lowheight = highheight
554	lowug_vert = highug_vert
555	lowvg_vert = highvg_vert
556	lowwsubs_vert = highwsubs_vert
557	low_td_lsa_lpt = high_td_lsa_lpt
558	low_td_lsa_q = high_td_lsa_q
559	low_td_sub_lpt = high_td_sub_lpt
560	low_td_sub_q = high_td_sub_q
561
562	ierrn = 0
563	READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert, &
564	highvg_vert, highwsubs_vert, &
565	high_td_lsa_lpt, &
566	high_td_lsa_q, &
567	high_td_sub_lpt, high_td_sub_q
568
569	IF ( ierrn /= 0 ) THEN
570	WRITE( message_string, * ) 'zu(',k,') = ', zu(k), 'm ', &
571	'is higher than the maximum height in LSF_DATA which ',&
572	'is ', lowheight, 'm. Interpolation on PALM ', &
573	'grid is not possible.'
574	CALL message( 'ls_forcing', 'PA0395', 1, 2, 0, 6, 0 )
575	ENDIF
576
577	ENDIF
578
579	!
580	!-- Interpolation of prescribed profiles in space
581	fac = (highheight-zu(k))/(highheight - lowheight)
582
583	ug_vert(k,nt) = fac * lowug_vert &
584	+ ( 1.0_wp - fac ) * highug_vert
585	vg_vert(k,nt) = fac * lowvg_vert &
586	+ ( 1.0_wp - fac ) * highvg_vert
587	wsubs_vert(k,nt) = fac * lowwsubs_vert &
588	+ ( 1.0_wp - fac ) * highwsubs_vert
589
590	td_lsa_lpt(k,nt) = fac * low_td_lsa_lpt &
591	+ ( 1.0_wp - fac ) * high_td_lsa_lpt
592	td_lsa_q(k,nt) = fac * low_td_lsa_q &
593	+ ( 1.0_wp - fac ) * high_td_lsa_q
594	td_sub_lpt(k,nt) = fac * low_td_sub_lpt &
595	+ ( 1.0_wp - fac ) * high_td_sub_lpt
596	td_sub_q(k,nt) = fac * low_td_sub_q &
597	+ ( 1.0_wp - fac ) * high_td_sub_q
598
599	ENDDO
600
601	ENDDO
602
603	!
604	!-- Large scale vertical velocity has to be zero at the surface
605	wsubs_vert(nzb,:) = 0.0_wp
606
607	IF ( time_vert(1) > end_time ) THEN
608	WRITE ( message_string, * ) 'Time dependent large scale profile ', &
609	'forcing from&LSF_DATA sets in after end of ' , &
610	'simulation - lsf_vert is set to FALSE'
611	CALL message( 'ls_forcing', 'PA0373', 0, 0, 0, 6, 0 )
612	lsf_vert = .FALSE.
613	ENDIF
614
615	CLOSE( finput )
616
617
618	END SUBROUTINE lsf_init
619
620	!------------------------------------------------------------------------------!
621	! Description:
622	! ------------
623	!> @todo Missing subroutine description.
624	!------------------------------------------------------------------------------!
625	SUBROUTINE ls_forcing_surf ( time )
626
627	IMPLICIT NONE
628
629	INTEGER(iwp) :: nt !<
630
631	REAL(wp) :: dum_surf_flux !<
632	REAL(wp) :: fac !<
633	REAL(wp), INTENT(in) :: time !<
634
635	!
636	!-- Interpolation in time of LSF_DATA at the surface
637	nt = 1
638	DO WHILE ( time > time_surf(nt) )
639	nt = nt + 1
640	ENDDO
641	IF ( time /= time_surf(nt) ) THEN
642	nt = nt - 1
643	ENDIF
644
645	fac = ( time -time_surf(nt) ) / ( time_surf(nt+1) - time_surf(nt) )
646
647	IF ( ibc_pt_b == 0 ) THEN
648	!
649	!-- In case of Dirichlet boundary condition shf must not
650	!-- be set - it is calculated via MOST in prandtl_fluxes
651	pt_surface = pt_surf(nt) + fac * ( pt_surf(nt+1) - pt_surf(nt) )
652
653	ELSEIF ( ibc_pt_b == 1 ) THEN
654	!
655	!-- In case of Neumann boundary condition pt_surface is needed for
656	!-- calculation of reference density
657	dum_surf_flux = ( shf_surf(nt) + fac * &
658	( shf_surf(nt+1) - shf_surf(nt) ) &
659	) * heatflux_input_conversion(nzb)
660	!
661	!-- Save surface sensible heat flux on default, natural and urban surface
662	!-- type, if required
663	IF ( surf_def_h(0)%ns >= 1 ) surf_def_h(0)%shf(:) = dum_surf_flux
664	IF ( surf_lsm_h%ns >= 1 ) surf_lsm_h%shf(:) = dum_surf_flux
665	IF ( surf_usm_h%ns >= 1 ) surf_usm_h%shf(:) = dum_surf_flux
666
667	pt_surface = pt_surf(nt) + fac * ( pt_surf(nt+1) - pt_surf(nt) )
668
669	ENDIF
670
671	IF ( ibc_q_b == 0 ) THEN
672	!
673	!-- In case of Dirichlet boundary condition qsws must not
674	!-- be set - it is calculated via MOST in prandtl_fluxes
675	q_surface = q_surf(nt) + fac * ( q_surf(nt+1) - q_surf(nt) )
676
677	ELSEIF ( ibc_q_b == 1 ) THEN
678	dum_surf_flux = ( qsws_surf(nt) + fac * &
679	( qsws_surf(nt+1) - qsws_surf(nt) ) &
680	) * waterflux_input_conversion(nzb)
681	!
682	!-- Save surface sensible heat flux on default, natural and urban surface
683	!-- type, if required
684	IF ( surf_def_h(0)%ns >= 1 ) surf_def_h(0)%qsws(:) = dum_surf_flux
685	IF ( surf_lsm_h%ns >= 1 ) surf_lsm_h%qsws(:) = dum_surf_flux
686	IF ( surf_usm_h%ns >= 1 ) surf_usm_h%qsws(:) = dum_surf_flux
687
688	ENDIF
689
690	surface_pressure = p_surf(nt) + fac * ( p_surf(nt+1) - p_surf(nt) )
691
692	END SUBROUTINE ls_forcing_surf
693
694
695
696
697	!------------------------------------------------------------------------------!
698	! Description:
699	! ------------
700	!> @todo Missing subroutine description.
701	!------------------------------------------------------------------------------!
702	SUBROUTINE ls_forcing_vert ( time )
703
704
705	IMPLICIT NONE
706
707	INTEGER(iwp) :: nt !<
708
709	REAL(wp) :: fac !<
710	REAL(wp), INTENT(in) :: time !<
711
712	!
713	!-- Interpolation in time of LSF_DATA for ug, vg and w_subs
714	nt = 1
715	DO WHILE ( time > time_vert(nt) )
716	nt = nt + 1
717	ENDDO
718	IF ( time /= time_vert(nt) ) THEN
719	nt = nt - 1
720	ENDIF
721
722	fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
723
724	ug = ug_vert(:,nt) + fac * ( ug_vert(:,nt+1) - ug_vert(:,nt) )
725	vg = vg_vert(:,nt) + fac * ( vg_vert(:,nt+1) - vg_vert(:,nt) )
726
727	IF ( large_scale_subsidence ) THEN
728	w_subs = wsubs_vert(:,nt) &
729	+ fac * ( wsubs_vert(:,nt+1) - wsubs_vert(:,nt) )
730	ENDIF
731
732	END SUBROUTINE ls_forcing_vert
733
734
735	!------------------------------------------------------------------------------!
736	! Description:
737	! ------------
738	!> Call for all grid points
739	!------------------------------------------------------------------------------!
740	SUBROUTINE ls_advec ( time, prog_var )
741
742
743	IMPLICIT NONE
744
745	CHARACTER (LEN=*) :: prog_var !<
746
747	REAL(wp), INTENT(in) :: time !<
748	REAL(wp) :: fac !<
749
750	INTEGER(iwp) :: i !<
751	INTEGER(iwp) :: j !<
752	INTEGER(iwp) :: k !<
753	INTEGER(iwp) :: nt !<
754
755	!
756	!-- Interpolation in time of LSF_DATA
757	nt = 1
758	DO WHILE ( time > time_vert(nt) )
759	nt = nt + 1
760	ENDDO
761	IF ( time /= time_vert(nt) ) THEN
762	nt = nt - 1
763	ENDIF
764
765	fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
766
767	!
768	!-- Add horizontal large scale advection tendencies of pt and q
769	SELECT CASE ( prog_var )
770
771	CASE ( 'pt' )
772
773	DO i = nxl, nxr
774	DO j = nys, nyn
775	DO k = nzb+1, nzt
776	tend(k,j,i) = tend(k,j,i) + td_lsa_lpt(k,nt) + fac * &
777	( td_lsa_lpt(k,nt+1) - td_lsa_lpt(k,nt) ) *&
778	MERGE( 1.0_wp, 0.0_wp, &
779	BTEST( wall_flags_0(k,j,i), 0 ) )
780	ENDDO
781	ENDDO
782	ENDDO
783
784	CASE ( 'q' )
785
786	DO i = nxl, nxr
787	DO j = nys, nyn
788	DO k = nzb+1, nzt
789	tend(k,j,i) = tend(k,j,i) + td_lsa_q(k,nt) + fac * &
790	( td_lsa_q(k,nt+1) - td_lsa_q(k,nt) ) * &
791	MERGE( 1.0_wp, 0.0_wp, &
792	BTEST( wall_flags_0(k,j,i), 0 ) )
793	ENDDO
794	ENDDO
795	ENDDO
796
797	END SELECT
798
799	!
800	!-- Subsidence of pt and q with prescribed subsidence tendencies
801	IF ( large_scale_subsidence .AND. use_subsidence_tendencies ) THEN
802
803	SELECT CASE ( prog_var )
804
805	CASE ( 'pt' )
806
807	DO i = nxl, nxr
808	DO j = nys, nyn
809	DO k = nzb+1, nzt
810	tend(k,j,i) = tend(k,j,i) + td_sub_lpt(k,nt) + fac * &
811	( td_sub_lpt(k,nt+1) - td_sub_lpt(k,nt) )*&
812	MERGE( 1.0_wp, 0.0_wp, &
813	BTEST( wall_flags_0(k,j,i), 0 ) )
814	ENDDO
815	ENDDO
816	ENDDO
817
818	CASE ( 'q' )
819
820	DO i = nxl, nxr
821	DO j = nys, nyn
822	DO k = nzb+1, nzt
823	tend(k,j,i) = tend(k,j,i) + td_sub_q(k,nt) + fac * &
824	( td_sub_q(k,nt+1) - td_sub_q(k,nt) ) * &
825	MERGE( 1.0_wp, 0.0_wp, &
826	BTEST( wall_flags_0(k,j,i), 0 ) )
827	ENDDO
828	ENDDO
829	ENDDO
830
831	END SELECT
832
833	ENDIF
834
835	END SUBROUTINE ls_advec
836
837
838	!------------------------------------------------------------------------------!
839	! Description:
840	! ------------
841	!> Call for grid point i,j
842	!------------------------------------------------------------------------------!
843	SUBROUTINE ls_advec_ij ( i, j, time, prog_var )
844
845	IMPLICIT NONE
846
847	CHARACTER (LEN=*) :: prog_var !<
848
849	REAL(wp), INTENT(in) :: time !<
850	REAL(wp) :: fac !<
851
852	INTEGER(iwp) :: i !<
853	INTEGER(iwp) :: j !<
854	INTEGER(iwp) :: k !<
855	INTEGER(iwp) :: nt !<
856
857	!
858	!-- Interpolation in time of LSF_DATA
859	nt = 1
860	DO WHILE ( time > time_vert(nt) )
861	nt = nt + 1
862	ENDDO
863	IF ( time /= time_vert(nt) ) THEN
864	nt = nt - 1
865	ENDIF
866
867	fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
868
869	!
870	!-- Add horizontal large scale advection tendencies of pt and q
871	SELECT CASE ( prog_var )
872
873	CASE ( 'pt' )
874
875	DO k = nzb+1, nzt
876	tend(k,j,i) = tend(k,j,i) + td_lsa_lpt(k,nt) &
877	+ fac * ( td_lsa_lpt(k,nt+1) - td_lsa_lpt(k,nt) )*&
878	MERGE( 1.0_wp, 0.0_wp, &
879	BTEST( wall_flags_0(k,j,i), 0 ) )
880	ENDDO
881
882	CASE ( 'q' )
883
884	DO k = nzb+1, nzt
885	tend(k,j,i) = tend(k,j,i) + td_lsa_q(k,nt) &
886	+ fac * ( td_lsa_q(k,nt+1) - td_lsa_q(k,nt) ) * &
887	MERGE( 1.0_wp, 0.0_wp, &
888	BTEST( wall_flags_0(k,j,i), 0 ) )
889	ENDDO
890
891	END SELECT
892
893	!
894	!-- Subsidence of pt and q with prescribed profiles
895	IF ( large_scale_subsidence .AND. use_subsidence_tendencies ) THEN
896
897	SELECT CASE ( prog_var )
898
899	CASE ( 'pt' )
900
901	DO k = nzb+1, nzt
902	tend(k,j,i) = tend(k,j,i) + td_sub_lpt(k,nt) + fac * &
903	( td_sub_lpt(k,nt+1) - td_sub_lpt(k,nt) ) * &
904	MERGE( 1.0_wp, 0.0_wp, &
905	BTEST( wall_flags_0(k,j,i), 0 ) )
906	ENDDO
907
908	CASE ( 'q' )
909
910	DO k = nzb+1, nzt
911	tend(k,j,i) = tend(k,j,i) + td_sub_q(k,nt) + fac * &
912	( td_sub_q(k,nt+1) - td_sub_q(k,nt) ) * &
913	MERGE( 1.0_wp, 0.0_wp, &
914	BTEST( wall_flags_0(k,j,i), 0 ) )
915	ENDDO
916
917	END SELECT
918
919	ENDIF
920
921	END SUBROUTINE ls_advec_ij
922
923
924	!------------------------------------------------------------------------------!
925	! Description:
926	! ------------
927	!> @todo Missing subroutine description.
928	!------------------------------------------------------------------------------!
929	SUBROUTINE nudge_init
930
931	IMPLICIT NONE
932
933
934	INTEGER(iwp) :: finput = 90 !<
935	INTEGER(iwp) :: ierrn !<
936	INTEGER(iwp) :: k !<
937	INTEGER(iwp) :: nt !<
938
939	CHARACTER(1) :: hash !<
940
941	REAL(wp) :: highheight !<
942	REAL(wp) :: highqnudge !<
943	REAL(wp) :: highptnudge !<
944	REAL(wp) :: highunudge !<
945	REAL(wp) :: highvnudge !<
946	REAL(wp) :: highwnudge !<
947	REAL(wp) :: hightnudge !<
948
949	REAL(wp) :: lowheight !<
950	REAL(wp) :: lowqnudge !<
951	REAL(wp) :: lowptnudge !<
952	REAL(wp) :: lowunudge !<
953	REAL(wp) :: lowvnudge !<
954	REAL(wp) :: lowwnudge !<
955	REAL(wp) :: lowtnudge !<
956
957	REAL(wp) :: fac !<
958
959	ALLOCATE( ptnudge(nzb:nzt+1,1:ntnudge), qnudge(nzb:nzt+1,1:ntnudge), &
960	tnudge(nzb:nzt+1,1:ntnudge), unudge(nzb:nzt+1,1:ntnudge), &
961	vnudge(nzb:nzt+1,1:ntnudge), wnudge(nzb:nzt+1,1:ntnudge) )
962
963	ALLOCATE( tmp_tnudge(nzb:nzt) )
964
965	ALLOCATE( timenudge(0:ntnudge) )
966
967	ptnudge = 0.0_wp; qnudge = 0.0_wp; tnudge = 0.0_wp; unudge = 0.0_wp
968	vnudge = 0.0_wp; wnudge = 0.0_wp; timenudge = 0.0_wp
969	!
970	!-- Initialize array tmp_nudge with a current nudging time scale of 6 hours
971	tmp_tnudge = 21600.0_wp
972
973	nt = 0
974	OPEN ( finput, FILE='NUDGING_DATA', STATUS='OLD', &
975	FORM='FORMATTED', IOSTAT=ierrn )
976
977	IF ( ierrn /= 0 ) THEN
978	message_string = 'file NUDGING_DATA does not exist'
979	CALL message( 'nudging', 'PA0365', 1, 2, 0, 6, 0 )
980	ENDIF
981
982	ierrn = 0
983
984	rloop:DO
985	nt = nt + 1
986	hash = "#"
987	ierrn = 1 ! not zero
988	!
989	!-- Search for the next line consisting of "# time",
990	!-- from there onwards the profiles will be read
991	DO WHILE ( .NOT. ( hash == "#" .AND. ierrn == 0 ) )
992
993	READ ( finput, *, IOSTAT=ierrn ) hash, timenudge(nt)
994	IF ( ierrn < 0 ) EXIT rloop
995
996	ENDDO
997
998	ierrn = 0
999	READ ( finput, *, IOSTAT=ierrn ) lowheight, lowtnudge, lowunudge, &
1000	lowvnudge, lowwnudge , lowptnudge, &
1001	lowqnudge
1002
1003	IF ( ierrn /= 0 ) THEN
1004	message_string = 'errors in file NUDGING_DATA'
1005	CALL message( 'nudging', 'PA0366', 1, 2, 0, 6, 0 )
1006	ENDIF
1007
1008	ierrn = 0
1009	READ ( finput, *, IOSTAT=ierrn ) highheight, hightnudge, highunudge, &
1010	highvnudge, highwnudge , highptnudge, &
1011	highqnudge
1012
1013	IF ( ierrn /= 0 ) THEN
1014	message_string = 'errors in file NUDGING_DATA'
1015	CALL message( 'nudging', 'PA0366', 1, 2, 0, 6, 0 )
1016	ENDIF
1017
1018	DO k = nzb, nzt+1
1019	DO WHILE ( highheight < zu(k) )
1020	lowheight = highheight
1021	lowtnudge = hightnudge
1022	lowunudge = highunudge
1023	lowvnudge = highvnudge
1024	lowwnudge = highwnudge
1025	lowptnudge = highptnudge
1026	lowqnudge = highqnudge
1027
1028	ierrn = 0
1029	READ ( finput, *, IOSTAT=ierrn ) highheight , hightnudge , &
1030	highunudge , highvnudge , &
1031	highwnudge , highptnudge, &
1032	highqnudge
1033	IF (ierrn /= 0 ) THEN
1034	WRITE( message_string, * ) 'zu(',k,') = ', zu(k), 'm is ',&
1035	'higher than the maximum height in NUDING_DATA which ', &
1036	'is ', lowheight, 'm. Interpolation on PALM ', &
1037	'grid is not possible.'
1038	CALL message( 'nudging', 'PA0364', 1, 2, 0, 6, 0 )
1039	ENDIF
1040	ENDDO
1041
1042	!
1043	!-- Interpolation of prescribed profiles in space
1044
1045	fac = ( highheight - zu(k) ) / ( highheight - lowheight )
1046
1047	tnudge(k,nt) = fac * lowtnudge + ( 1.0_wp - fac ) * hightnudge
1048	unudge(k,nt) = fac * lowunudge + ( 1.0_wp - fac ) * highunudge
1049	vnudge(k,nt) = fac * lowvnudge + ( 1.0_wp - fac ) * highvnudge
1050	wnudge(k,nt) = fac * lowwnudge + ( 1.0_wp - fac ) * highwnudge
1051	ptnudge(k,nt) = fac * lowptnudge + ( 1.0_wp - fac ) * highptnudge
1052	qnudge(k,nt) = fac * lowqnudge + ( 1.0_wp - fac ) * highqnudge
1053	ENDDO
1054
1055	ENDDO rloop
1056
1057	CLOSE ( finput )
1058
1059	!
1060	!-- Overwrite initial profiles in case of nudging
1061	IF ( nudging ) THEN
1062	pt_init = ptnudge(:,1)
1063	u_init = unudge(:,1)
1064	v_init = vnudge(:,1)
1065	IF ( humidity ) THEN ! is passive_scalar correct???
1066	q_init = qnudge(:,1)
1067	ENDIF
1068
1069	WRITE( message_string, * ) 'Initial profiles of u, v and ', &
1070	'scalars from NUDGING_DATA are used.'
1071	CALL message( 'large_scale_forcing_nudging', 'PA0370', 0, 0, 0, 6, 0 )
1072	ENDIF
1073
1074
1075	END SUBROUTINE nudge_init
1076
1077	!------------------------------------------------------------------------------!
1078	! Description:
1079	! ------------
1080	!> @todo Missing subroutine description.
1081	!------------------------------------------------------------------------------!
1082	SUBROUTINE calc_tnudge ( time )
1083
1084	IMPLICIT NONE
1085
1086
1087	REAL(wp) :: dtm !<
1088	REAL(wp) :: dtp !<
1089	REAL(wp) :: time !<
1090
1091	INTEGER(iwp) :: k !<
1092	INTEGER(iwp) :: nt !<
1093
1094	nt = 1
1095	DO WHILE ( time > timenudge(nt) )
1096	nt = nt+1
1097	ENDDO
1098	IF ( time /= timenudge(1) ) THEN
1099	nt = nt-1
1100	ENDIF
1101
1102	dtm = ( time - timenudge(nt) ) / ( timenudge(nt+1) - timenudge(nt) )
1103	dtp = ( timenudge(nt+1) - time ) / ( timenudge(nt+1) - timenudge(nt) )
1104
1105	DO k = nzb, nzt
1106	tmp_tnudge(k) = MAX( dt_3d, tnudge(k,nt) * dtp + tnudge(k,nt+1) * dtm )
1107	ENDDO
1108
1109	END SUBROUTINE calc_tnudge
1110
1111	!------------------------------------------------------------------------------!
1112	! Description:
1113	! ------------
1114	!> Call for all grid points
1115	!------------------------------------------------------------------------------!
1116	SUBROUTINE nudge ( time, prog_var )
1117
1118	IMPLICIT NONE
1119
1120	CHARACTER (LEN=*) :: prog_var !<
1121
1122	REAL(wp) :: tmp_tend !<
1123	REAL(wp) :: dtm !<
1124	REAL(wp) :: dtp !<
1125	REAL(wp) :: time !<
1126
1127	INTEGER(iwp) :: i !<
1128	INTEGER(iwp) :: j !<
1129	INTEGER(iwp) :: k !<
1130	INTEGER(iwp) :: nt !<
1131
1132
1133	nt = 1
1134	DO WHILE ( time > timenudge(nt) )
1135	nt = nt+1
1136	ENDDO
1137	IF ( time /= timenudge(1) ) THEN
1138	nt = nt-1
1139	ENDIF
1140
1141	dtm = ( time - timenudge(nt) ) / ( timenudge(nt+1) - timenudge(nt) )
1142	dtp = ( timenudge(nt+1) - time ) / ( timenudge(nt+1) - timenudge(nt) )
1143
1144	SELECT CASE ( prog_var )
1145
1146	CASE ( 'u' )
1147
1148	DO i = nxl, nxr
1149	DO j = nys, nyn
1150
1151	DO k = nzb+1, nzt
1152
1153	tmp_tend = - ( hom(k,1,1,0) - ( unudge(k,nt) * dtp + &
1154	unudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
1155
1156	tend(k,j,i) = tend(k,j,i) + tmp_tend * &
1157	MERGE( 1.0_wp, 0.0_wp, &
1158	BTEST( wall_flags_0(k,j,i), 1 ) )
1159
1160	sums_ls_l(k,6) = sums_ls_l(k,6) + tmp_tend * &
1161	weight_substep(intermediate_timestep_count)
1162	ENDDO
1163
1164	sums_ls_l(nzt+1,6) = sums_ls_l(nzt,6)
1165
1166	ENDDO
1167	ENDDO
1168
1169	CASE ( 'v' )
1170
1171	DO i = nxl, nxr
1172	DO j = nys, nyn
1173
1174	DO k = nzb+1, nzt
1175
1176	tmp_tend = - ( hom(k,1,2,0) - ( vnudge(k,nt) * dtp + &
1177	vnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
1178
1179	tend(k,j,i) = tend(k,j,i) + tmp_tend * &
1180	MERGE( 1.0_wp, 0.0_wp, &
1181	BTEST( wall_flags_0(k,j,i), 2 ) )
1182
1183	sums_ls_l(k,7) = sums_ls_l(k,7) + tmp_tend * &
1184	weight_substep(intermediate_timestep_count)
1185	ENDDO
1186
1187	sums_ls_l(nzt+1,7) = sums_ls_l(nzt,7)
1188
1189	ENDDO
1190	ENDDO
1191
1192	CASE ( 'pt' )
1193
1194	DO i = nxl, nxr
1195	DO j = nys, nyn
1196
1197	DO k = nzb+1, nzt
1198
1199	tmp_tend = - ( hom(k,1,4,0) - ( ptnudge(k,nt) * dtp + &
1200	ptnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
1201
1202	tend(k,j,i) = tend(k,j,i) + tmp_tend * &
1203	MERGE( 1.0_wp, 0.0_wp, &
1204	BTEST( wall_flags_0(k,j,i), 0 ) )
1205
1206	sums_ls_l(k,4) = sums_ls_l(k,4) + tmp_tend * &
1207	weight_substep(intermediate_timestep_count)
1208	ENDDO
1209
1210	sums_ls_l(nzt+1,4) = sums_ls_l(nzt,4)
1211
1212	ENDDO
1213	ENDDO
1214
1215	CASE ( 'q' )
1216
1217	DO i = nxl, nxr
1218	DO j = nys, nyn
1219
1220	DO k = nzb+1, nzt
1221
1222	tmp_tend = - ( hom(k,1,41,0) - ( qnudge(k,nt) * dtp + &
1223	qnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
1224
1225	tend(k,j,i) = tend(k,j,i) + tmp_tend * &
1226	MERGE( 1.0_wp, 0.0_wp, &
1227	BTEST( wall_flags_0(k,j,i), 0 ) )
1228
1229	sums_ls_l(k,5) = sums_ls_l(k,5) + tmp_tend * &
1230	weight_substep(intermediate_timestep_count)
1231	ENDDO
1232
1233	sums_ls_l(nzt+1,5) = sums_ls_l(nzt,5)
1234
1235	ENDDO
1236	ENDDO
1237
1238	CASE DEFAULT
1239	message_string = 'unknown prognostic variable "' // prog_var // '"'
1240	CALL message( 'nudge', 'PA0367', 1, 2, 0, 6, 0 )
1241
1242	END SELECT
1243
1244	END SUBROUTINE nudge
1245
1246
1247	!------------------------------------------------------------------------------!
1248	! Description:
1249	! ------------
1250	!> Call for grid point i,j
1251	!------------------------------------------------------------------------------!
1252
1253	SUBROUTINE nudge_ij( i, j, time, prog_var )
1254
1255	IMPLICIT NONE
1256
1257
1258	CHARACTER (LEN=*) :: prog_var !<
1259
1260	REAL(wp) :: tmp_tend !<
1261	REAL(wp) :: dtm !<
1262	REAL(wp) :: dtp !<
1263	REAL(wp) :: time !<
1264
1265	INTEGER(iwp) :: i !<
1266	INTEGER(iwp) :: j !<
1267	INTEGER(iwp) :: k !<
1268	INTEGER(iwp) :: nt !<
1269
1270
1271	nt = 1
1272	DO WHILE ( time > timenudge(nt) )
1273	nt = nt+1
1274	ENDDO
1275	IF ( time /= timenudge(1) ) THEN
1276	nt = nt-1
1277	ENDIF
1278
1279	dtm = ( time - timenudge(nt) ) / ( timenudge(nt+1) - timenudge(nt) )
1280	dtp = ( timenudge(nt+1) - time ) / ( timenudge(nt+1) - timenudge(nt) )
1281
1282	SELECT CASE ( prog_var )
1283
1284	CASE ( 'u' )
1285
1286	DO k = nzb+1, nzt
1287
1288	tmp_tend = - ( hom(k,1,1,0) - ( unudge(k,nt) * dtp + &
1289	unudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
1290
1291	tend(k,j,i) = tend(k,j,i) + tmp_tend * &
1292	MERGE( 1.0_wp, 0.0_wp, &
1293	BTEST( wall_flags_0(k,j,i), 1 ) )
1294
1295	sums_ls_l(k,6) = sums_ls_l(k,6) + tmp_tend &
1296	* weight_substep(intermediate_timestep_count)
1297	ENDDO
1298
1299	sums_ls_l(nzt+1,6) = sums_ls_l(nzt,6)
1300
1301	CASE ( 'v' )
1302
1303	DO k = nzb+1, nzt
1304
1305	tmp_tend = - ( hom(k,1,2,0) - ( vnudge(k,nt) * dtp + &
1306	vnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
1307
1308	tend(k,j,i) = tend(k,j,i) + tmp_tend * &
1309	MERGE( 1.0_wp, 0.0_wp, &
1310	BTEST( wall_flags_0(k,j,i), 2 ) )
1311
1312	sums_ls_l(k,7) = sums_ls_l(k,7) + tmp_tend &
1313	* weight_substep(intermediate_timestep_count)
1314	ENDDO
1315
1316	sums_ls_l(nzt+1,7) = sums_ls_l(nzt,7)
1317
1318	CASE ( 'pt' )
1319
1320	DO k = nzb+1, nzt
1321
1322	tmp_tend = - ( hom(k,1,4,0) - ( ptnudge(k,nt) * dtp + &
1323	ptnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
1324
1325	tend(k,j,i) = tend(k,j,i) + tmp_tend * &
1326	MERGE( 1.0_wp, 0.0_wp, &
1327	BTEST( wall_flags_0(k,j,i), 0 ) )
1328
1329	sums_ls_l(k,4) = sums_ls_l(k,4) + tmp_tend &
1330	* weight_substep(intermediate_timestep_count)
1331	ENDDO
1332
1333	sums_ls_l(nzt+1,4) = sums_ls_l(nzt,4)
1334
1335
1336	CASE ( 'q' )
1337
1338	DO k = nzb+1, nzt
1339
1340	tmp_tend = - ( hom(k,1,41,0) - ( qnudge(k,nt) * dtp + &
1341	qnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
1342
1343	tend(k,j,i) = tend(k,j,i) + tmp_tend * &
1344	MERGE( 1.0_wp, 0.0_wp, &
1345	BTEST( wall_flags_0(k,j,i), 0 ) )
1346
1347	sums_ls_l(k,5) = sums_ls_l(k,5) + tmp_tend &
1348	* weight_substep(intermediate_timestep_count)
1349	ENDDO
1350
1351	sums_ls_l(nzt+1,5) = sums_ls_l(nzt,5)
1352
1353	CASE DEFAULT
1354	message_string = 'unknown prognostic variable "' // prog_var // '"'
1355	CALL message( 'nudge', 'PA0367', 1, 2, 0, 6, 0 )
1356
1357	END SELECT
1358
1359
1360	END SUBROUTINE nudge_ij
1361
1362
1363	!------------------------------------------------------------------------------!
1364	! Description:
1365	! ------------
1366	!> @todo Missing subroutine description.
1367	!------------------------------------------------------------------------------!
1368	SUBROUTINE nudge_ref ( time )
1369
1370	IMPLICIT NONE
1371
1372	INTEGER(iwp) :: nt !<
1373
1374	REAL(wp) :: fac !<
1375	REAL(wp), INTENT(in) :: time !<
1376
1377	!
1378	!-- Interpolation in time of NUDGING_DATA for pt_init and q_init. This is
1379	!-- needed for correct upper boundary conditions for pt and q and in case that
1380	! large scale subsidence as well as scalar Rayleigh-damping are used
1381	nt = 1
1382	DO WHILE ( time > time_vert(nt) )
1383	nt = nt + 1
1384	ENDDO
1385	IF ( time /= time_vert(nt) ) THEN
1386	nt = nt - 1
1387	ENDIF
1388
1389	fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
1390
1391	pt_init = ptnudge(:,nt) + fac * ( ptnudge(:,nt+1) - ptnudge(:,nt) )
1392	q_init = qnudge(:,nt) + fac * ( qnudge(:,nt+1) - qnudge(:,nt) )
1393	u_init = unudge(:,nt) + fac * ( unudge(:,nt+1) - unudge(:,nt) )
1394	v_init = vnudge(:,nt) + fac * ( vnudge(:,nt+1) - vnudge(:,nt) )
1395
1396	END SUBROUTINE nudge_ref
1397
1398	END MODULE lsf_nudging_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |