Home

Context Navigation

large_scale_forcing_nudging_mod.f90 @ 4858

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

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

Context Navigation

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

Download in other formats: