Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |