Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |