Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |