Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |