Home

Context Navigation

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

Last change on this file since 3424 was 3347, checked in by suehring, 6 years ago
Offline nesting revised and separated from large_scale_forcing_mod; Time-dependent synthetic turbulence generator; bugfixes in USM/LSM radiation model and init_pegrid
Property svn:keywords set to `Id`
File size: 55.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |