Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |