Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |