Home

Context Navigation

source: palm/trunk/SOURCE/model_1d_mod.f90 @ 2943

Last change on this file since 2943 was 2918, checked in by gronemeier, 7 years ago
moved initialization of mixing length to turbulence_closure_mod; consider walls in calculation of ml in rans-mode
Property svn:keywords set to `Id`
File size: 42.1 KB

Line
1	!> @file model_1d_mod.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of the PALM model system.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the
6	! terms of the GNU General Public License as published by the Free Software
7	! Foundation, either version 3 of the License, or (at your option) any later
8	! version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 1997-2018 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: model_1d_mod.f90 2918 2018-03-21 15:52:14Z suehring $
27	! - rename l_black into l1d_init
28	! - calculate l_grid within init_1d_model and save it as l1d_init
29	!
30	! 2718 2018-01-02 08:49:38Z maronga
31	! Corrected "Former revisions" section
32	!
33	! 2696 2017-12-14 17:12:51Z kanani
34	! Change in file header (GPL part)
35	! implement TKE-e closure
36	! modification of dissipation production according to Detering and Etling
37	! reduced factor for timestep criterion to 0.125 and first dt to 1s (TG)
38	!
39	! 2339 2017-08-07 13:55:26Z gronemeier
40	! corrected timestamp in header
41	!
42	! 2338 2017-08-07 12:15:38Z gronemeier
43	! renamed init_1d_model to model_1d_mod and and formatted it as a module;
44	! reformatted output of profiles
45	!
46	! 2337 2017-08-07 08:59:53Z gronemeier
47	! revised calculation of mixing length
48	! removed rounding of time step
49	! corrected calculation of virtual potential temperature
50	!
51	! 2334 2017-08-04 11:57:04Z gronemeier
52	! set c_m = 0.4 according to Detering and Etling (1985)
53	!
54	! 2299 2017-06-29 10:14:38Z maronga
55	! Removed german text
56	!
57	! 2101 2017-01-05 16:42:31Z suehring
58	!
59	! 2059 2016-11-10 14:20:40Z maronga
60	! Corrected min/max values of Rif.
61	!
62	! 2000 2016-08-20 18:09:15Z knoop
63	! Forced header and separation lines into 80 columns
64	!
65	! 1960 2016-07-12 16:34:24Z suehring
66	! Remove passive_scalar from IF-statements, as 1D-scalar profile is effectively
67	! not used.
68	! Formatting adjustment
69	!
70	! 1808 2016-04-05 19:44:00Z raasch
71	! routine local_flush replaced by FORTRAN statement
72	!
73	! 1709 2015-11-04 14:47:01Z maronga
74	! Set initial time step to 10 s to avoid instability of the 1d model for small
75	! grid spacings
76	!
77	! 1697 2015-10-28 17:14:10Z raasch
78	! small E- and F-FORMAT changes to avoid informative compiler messages about
79	! insufficient field width
80	!
81	! 1691 2015-10-26 16:17:44Z maronga
82	! Renamed prandtl_layer to constant_flux_layer. rif is replaced by ol and zeta.
83	!
84	! 1682 2015-10-07 23:56:08Z knoop
85	! Code annotations made doxygen readable
86	!
87	! 1353 2014-04-08 15:21:23Z heinze
88	! REAL constants provided with KIND-attribute
89	!
90	! 1346 2014-03-27 13:18:20Z heinze
91	! Bugfix: REAL constants provided with KIND-attribute especially in call of
92	! intrinsic function like MAX, MIN, SIGN
93	!
94	! 1322 2014-03-20 16:38:49Z raasch
95	! REAL functions provided with KIND-attribute
96	!
97	! 1320 2014-03-20 08:40:49Z raasch
98	! ONLY-attribute added to USE-statements,
99	! kind-parameters added to all INTEGER and REAL declaration statements,
100	! kinds are defined in new module kinds,
101	! revision history before 2012 removed,
102	! comment fields (!:) to be used for variable explanations added to
103	! all variable declaration statements
104	!
105	! 1036 2012-10-22 13:43:42Z raasch
106	! code put under GPL (PALM 3.9)
107	!
108	! 1015 2012-09-27 09:23:24Z raasch
109	! adjustment of mixing length to the Prandtl mixing length at first grid point
110	! above ground removed
111	!
112	! 1001 2012-09-13 14:08:46Z raasch
113	! all actions concerning leapfrog scheme removed
114	!
115	! 996 2012-09-07 10:41:47Z raasch
116	! little reformatting
117	!
118	! 978 2012-08-09 08:28:32Z fricke
119	! roughness length for scalar quantities z0h1d added
120	!
121	! Revision 1.1 1998/03/09 16:22:10 raasch
122	! Initial revision
123	!
124	!
125	! Description:
126	! ------------
127	!> 1D-model to initialize the 3D-arrays.
128	!> The temperature profile is set as steady and a corresponding steady solution
129	!> of the wind profile is being computed.
130	!> All subroutines required can be found within this file.
131	!>
132	!> @todo harmonize code with new surface_layer_fluxes module
133	!> @bug 1D model crashes when using small grid spacings in the order of 1 m
134	!> @fixme option "as_in_3d_model" seems to eb an inappropriate option because
135	!> the 1D model uses different turbulence closure approaches at least if
136	!> the 3D model is set to LES-mode.
137	!------------------------------------------------------------------------------!
138	MODULE model_1d_mod
139
140	USE arrays_3d, &
141	ONLY: dd2zu, ddzu, ddzw, dzu, dzw, pt_init, q_init, ug, u_init, &
142	vg, v_init, zu
143
144	USE control_parameters, &
145	ONLY: constant_diffusion, constant_flux_layer, dissipation_1d, f, g, &
146	humidity, ibc_e_b, intermediate_timestep_count, &
147	intermediate_timestep_count_max, kappa, km_constant, &
148	message_string, mixing_length_1d, prandtl_number, &
149	roughness_length, run_description_header, simulated_time_chr, &
150	timestep_scheme, tsc, z0h_factor
151
152	USE indices, &
153	ONLY: nzb, nzb_diff, nzt
154
155	USE kinds
156
157	USE pegrid
158
159
160	IMPLICIT NONE
161
162	INTEGER(iwp) :: current_timestep_number_1d = 0 !< current timestep number (1d-model)
163	INTEGER(iwp) :: damp_level_ind_1d !< lower grid index of damping layer (1d-model)
164
165	LOGICAL :: run_control_header_1d = .FALSE. !< flag for output of run control header (1d-model)
166	LOGICAL :: stop_dt_1d = .FALSE. !< termination flag, used in case of too small timestep (1d-model)
167
168	REAL(wp) :: c_1 = 1.44_wp !< model constant
169	REAL(wp) :: c_2 = 1.92_wp !< model constant
170	REAL(wp) :: c_3 = 1.44_wp !< model constant
171	REAL(wp) :: c_h = 0.0015_wp !< model constant according to Detering and Etling (1985)
172	REAL(wp) :: c_m = 0.4_wp !< model constant, 0.4 according to Detering and Etling (1985)
173	REAL(wp) :: c_mu = 0.09_wp !< model constant
174	REAL(wp) :: damp_level_1d = -1.0_wp !< namelist parameter
175	REAL(wp) :: dt_1d = 60.0_wp !< dynamic timestep (1d-model)
176	REAL(wp) :: dt_max_1d = 300.0_wp !< timestep limit (1d-model)
177	REAL(wp) :: dt_pr_1d = 9999999.9_wp !< namelist parameter
178	REAL(wp) :: dt_run_control_1d = 60.0_wp !< namelist parameter
179	REAL(wp) :: end_time_1d = 864000.0_wp !< namelist parameter
180	REAL(wp) :: qs1d !< characteristic humidity scale (1d-model)
181	REAL(wp) :: simulated_time_1d = 0.0_wp !< updated simulated time (1d-model)
182	REAL(wp) :: sig_diss = 1.3_wp !< model constant
183	REAL(wp) :: time_pr_1d = 0.0_wp !< updated simulated time for profile output (1d-model)
184	REAL(wp) :: time_run_control_1d = 0.0_wp !< updated simulated time for run-control output (1d-model)
185	REAL(wp) :: ts1d !< characteristic temperature scale (1d-model)
186	REAL(wp) :: us1d !< friction velocity (1d-model)
187	REAL(wp) :: usws1d !< u-component of the momentum flux (1d-model)
188	REAL(wp) :: vsws1d !< v-component of the momentum flux (1d-model)
189	REAL(wp) :: z01d !< roughness length for momentum (1d-model)
190	REAL(wp) :: z0h1d !< roughness length for scalars (1d-model)
191
192
193	REAL(wp), DIMENSION(:), ALLOCATABLE :: diss1d !< tke dissipation rate (1d-model)
194	REAL(wp), DIMENSION(:), ALLOCATABLE :: diss1d_p !< prognostic value of tke dissipation rate (1d-model)
195	REAL(wp), DIMENSION(:), ALLOCATABLE :: e1d !< tke (1d-model)
196	REAL(wp), DIMENSION(:), ALLOCATABLE :: e1d_p !< prognostic value of tke (1d-model)
197	REAL(wp), DIMENSION(:), ALLOCATABLE :: kh1d !< turbulent diffusion coefficient for heat (1d-model)
198	REAL(wp), DIMENSION(:), ALLOCATABLE :: km1d !< turbulent diffusion coefficient for momentum (1d-model)
199	REAL(wp), DIMENSION(:), ALLOCATABLE :: l1d !< mixing length for turbulent diffusion coefficients (1d-model)
200	REAL(wp), DIMENSION(:), ALLOCATABLE :: l1d_init !< initial mixing length (1d-model)
201	REAL(wp), DIMENSION(:), ALLOCATABLE :: l1d_diss !< mixing length for dissipation (1d-model)
202	REAL(wp), DIMENSION(:), ALLOCATABLE :: rif1d !< Richardson flux number (1d-model)
203	REAL(wp), DIMENSION(:), ALLOCATABLE :: te_diss !< tendency of diss (1d-model)
204	REAL(wp), DIMENSION(:), ALLOCATABLE :: te_dissm !< weighted tendency of diss for previous sub-timestep (1d-model)
205	REAL(wp), DIMENSION(:), ALLOCATABLE :: te_e !< tendency of e (1d-model)
206	REAL(wp), DIMENSION(:), ALLOCATABLE :: te_em !< weighted tendency of e for previous sub-timestep (1d-model)
207	REAL(wp), DIMENSION(:), ALLOCATABLE :: te_u !< tendency of u (1d-model)
208	REAL(wp), DIMENSION(:), ALLOCATABLE :: te_um !< weighted tendency of u for previous sub-timestep (1d-model)
209	REAL(wp), DIMENSION(:), ALLOCATABLE :: te_v !< tendency of v (1d-model)
210	REAL(wp), DIMENSION(:), ALLOCATABLE :: te_vm !< weighted tendency of v for previous sub-timestep (1d-model)
211	REAL(wp), DIMENSION(:), ALLOCATABLE :: u1d !< u-velocity component (1d-model)
212	REAL(wp), DIMENSION(:), ALLOCATABLE :: u1d_p !< prognostic value of u-velocity component (1d-model)
213	REAL(wp), DIMENSION(:), ALLOCATABLE :: v1d !< v-velocity component (1d-model)
214	REAL(wp), DIMENSION(:), ALLOCATABLE :: v1d_p !< prognostic value of v-velocity component (1d-model)
215
216	!
217	!-- Initialize 1D model
218	INTERFACE init_1d_model
219	MODULE PROCEDURE init_1d_model
220	END INTERFACE init_1d_model
221
222	!
223	!-- Print profiles
224	INTERFACE print_1d_model
225	MODULE PROCEDURE print_1d_model
226	END INTERFACE print_1d_model
227
228	!
229	!-- Print run control information
230	INTERFACE run_control_1d
231	MODULE PROCEDURE run_control_1d
232	END INTERFACE run_control_1d
233
234	!
235	!-- Main procedure
236	INTERFACE time_integration_1d
237	MODULE PROCEDURE time_integration_1d
238	END INTERFACE time_integration_1d
239
240	!
241	!-- Calculate time step
242	INTERFACE timestep_1d
243	MODULE PROCEDURE timestep_1d
244	END INTERFACE timestep_1d
245
246	SAVE
247
248	PRIVATE
249	!
250	!-- Public interfaces
251	PUBLIC init_1d_model
252
253	!
254	!-- Public variables
255	PUBLIC damp_level_1d, damp_level_ind_1d, diss1d, dt_pr_1d, &
256	dt_run_control_1d, e1d, end_time_1d, kh1d, km1d, l1d, rif1d, u1d, &
257	us1d, usws1d, v1d, vsws1d
258
259
260	CONTAINS
261
262	SUBROUTINE init_1d_model
263
264	USE grid_variables, &
265	ONLY: dx, dy
266
267	IMPLICIT NONE
268
269	CHARACTER (LEN=9) :: time_to_string !< function to transform time from real to character string
270
271	INTEGER(iwp) :: k !< loop index
272
273	REAL(wp) :: lambda !< maximum mixing length
274
275	!
276	!-- Allocate required 1D-arrays
277	ALLOCATE( diss1d(nzb:nzt+1), diss1d_p(nzb:nzt+1), &
278	e1d(nzb:nzt+1), e1d_p(nzb:nzt+1), kh1d(nzb:nzt+1), &
279	km1d(nzb:nzt+1), l1d(nzb:nzt+1), l1d_init(nzb:nzt+1), &
280	l1d_diss(nzb:nzt+1), rif1d(nzb:nzt+1), te_diss(nzb:nzt+1), &
281	te_dissm(nzb:nzt+1), te_e(nzb:nzt+1), &
282	te_em(nzb:nzt+1), te_u(nzb:nzt+1), te_um(nzb:nzt+1), &
283	te_v(nzb:nzt+1), te_vm(nzb:nzt+1), u1d(nzb:nzt+1), &
284	u1d_p(nzb:nzt+1), v1d(nzb:nzt+1), v1d_p(nzb:nzt+1) )
285
286	!
287	!-- Initialize arrays
288	IF ( constant_diffusion ) THEN
289	km1d = km_constant
290	kh1d = km_constant / prandtl_number
291	ELSE
292	diss1d = 0.0_wp; diss1d_p = 0.0_wp
293	e1d = 0.0_wp; e1d_p = 0.0_wp
294	kh1d = 0.0_wp; km1d = 0.0_wp
295	rif1d = 0.0_wp
296	!
297	!-- Compute the mixing length
298	l1d_init(nzb) = 0.0_wp
299
300	IF ( TRIM( mixing_length_1d ) == 'blackadar' ) THEN
301	!
302	!-- Blackadar mixing length
303	IF ( f /= 0.0_wp ) THEN
304	lambda = 2.7E-4_wp * SQRT( ug(nzt+1)2 + vg(nzt+1)2 ) / &
305	ABS( f ) + 1E-10_wp
306	ELSE
307	lambda = 30.0_wp
308	ENDIF
309
310	DO k = nzb+1, nzt+1
311	l1d_init(k) = kappa * zu(k) / ( 1.0_wp + kappa * zu(k) / lambda )
312	ENDDO
313
314	ELSEIF ( TRIM( mixing_length_1d ) == 'as_in_3d_model' ) THEN
315	!
316	!-- Use the same mixing length as in 3D model (LES-mode)
317	!### TODO: rename (delete?) this option
318	!### As the mixing length is different between RANS and LES mode, it
319	!### must be distinguished here between these modes. For RANS mode,
320	!### the mixing length is calculated accoding to Blackadar, which is
321	!### the other option at this point.
322	!### Maybe delete this option entirely (not appropriate in LES case)
323	!### 2018-03-20, gronemeier
324	DO k = nzb+1, nzt
325	l1d_init(k) = ( dx * dy * dzw(k) )**0.33333333333333_wp
326	ENDDO
327	l1d_init(nzt+1) = l1d_init(nzt)
328
329	ENDIF
330	ENDIF
331	l1d = l1d_init
332	l1d_diss = l1d_init
333	u1d = u_init
334	u1d_p = u_init
335	v1d = v_init
336	v1d_p = v_init
337
338	!
339	!-- Set initial horizontal velocities at the lowest grid levels to a very small
340	!-- value in order to avoid too small time steps caused by the diffusion limit
341	!-- in the initial phase of a run (at k=1, dz/2 occurs in the limiting formula!)
342	u1d(0:1) = 0.1_wp
343	u1d_p(0:1) = 0.1_wp
344	v1d(0:1) = 0.1_wp
345	v1d_p(0:1) = 0.1_wp
346
347	!
348	!-- For u, theta and the momentum fluxes plausible values are set
349	IF ( constant_flux_layer ) THEN
350	us1d = 0.1_wp ! without initial friction the flow would not change
351	ELSE
352	diss1d(nzb+1) = 1.0_wp
353	e1d(nzb+1) = 1.0_wp
354	km1d(nzb+1) = 1.0_wp
355	us1d = 0.0_wp
356	ENDIF
357	ts1d = 0.0_wp
358	usws1d = 0.0_wp
359	vsws1d = 0.0_wp
360	z01d = roughness_length
361	z0h1d = z0h_factor * z01d
362	IF ( humidity ) qs1d = 0.0_wp
363
364	!
365	!-- Tendencies must be preset in order to avoid runtime errors within the
366	!-- first Runge-Kutta step
367	te_dissm = 0.0_wp
368	te_em = 0.0_wp
369	te_um = 0.0_wp
370	te_vm = 0.0_wp
371
372	!
373	!-- Set model constant
374	IF ( dissipation_1d == 'as_in_3d_model' ) c_m = 0.1_wp
375
376	!
377	!-- Set start time in hh:mm:ss - format
378	simulated_time_chr = time_to_string( simulated_time_1d )
379
380	!
381	!-- Integrate the 1D-model equations using the Runge-Kutta scheme
382	CALL time_integration_1d
383
384
385	END SUBROUTINE init_1d_model
386
387
388
389	!------------------------------------------------------------------------------!
390	! Description:
391	! ------------
392	!> Runge-Kutta time differencing scheme for the 1D-model.
393	!------------------------------------------------------------------------------!
394
395	SUBROUTINE time_integration_1d
396
397	IMPLICIT NONE
398
399	CHARACTER (LEN=9) :: time_to_string !< function to transform time from real to character string
400
401	INTEGER(iwp) :: k !< loop index
402
403	REAL(wp) :: a !< auxiliary variable
404	REAL(wp) :: b !< auxiliary variable
405	REAL(wp) :: dpt_dz !< vertical temperature gradient
406	REAL(wp) :: flux !< vertical temperature gradient
407	REAL(wp) :: kmzm !< Km(z-dz/2)
408	REAL(wp) :: kmzp !< Km(z+dz/2)
409	REAL(wp) :: l_stable !< mixing length for stable case
410	REAL(wp) :: pt_0 !< reference temperature
411	REAL(wp) :: uv_total !< horizontal wind speed
412
413	!
414	!-- Determine the time step at the start of a 1D-simulation and
415	!-- determine and printout quantities used for run control
416	dt_1d = 1.0_wp
417	CALL run_control_1d
418
419	!
420	!-- Start of time loop
421	DO WHILE ( simulated_time_1d < end_time_1d .AND. .NOT. stop_dt_1d )
422
423	!
424	!-- Depending on the timestep scheme, carry out one or more intermediate
425	!-- timesteps
426
427	intermediate_timestep_count = 0
428	DO WHILE ( intermediate_timestep_count < &
429	intermediate_timestep_count_max )
430
431	intermediate_timestep_count = intermediate_timestep_count + 1
432
433	CALL timestep_scheme_steering
434
435	!
436	!-- Compute all tendency terms. If a constant-flux layer is simulated,
437	!-- k starts at nzb+2.
438	DO k = nzb_diff, nzt
439
440	kmzm = 0.5_wp * ( km1d(k-1) + km1d(k) )
441	kmzp = 0.5_wp * ( km1d(k) + km1d(k+1) )
442	!
443	!-- u-component
444	te_u(k) = f * ( v1d(k) - vg(k) ) + ( &
445	kmzp * ( u1d(k+1) - u1d(k) ) * ddzu(k+1) &
446	- kmzm * ( u1d(k) - u1d(k-1) ) * ddzu(k) &
447	) * ddzw(k)
448	!
449	!-- v-component
450	te_v(k) = -f * ( u1d(k) - ug(k) ) + ( &
451	kmzp * ( v1d(k+1) - v1d(k) ) * ddzu(k+1) &
452	- kmzm * ( v1d(k) - v1d(k-1) ) * ddzu(k) &
453	) * ddzw(k)
454	ENDDO
455	IF ( .NOT. constant_diffusion ) THEN
456	DO k = nzb_diff, nzt
457	!
458	!-- TKE and dissipation rate
459	kmzm = 0.5_wp * ( km1d(k-1) + km1d(k) )
460	kmzp = 0.5_wp * ( km1d(k) + km1d(k+1) )
461	IF ( .NOT. humidity ) THEN
462	pt_0 = pt_init(k)
463	flux = ( pt_init(k+1)-pt_init(k-1) ) * dd2zu(k)
464	ELSE
465	pt_0 = pt_init(k) * ( 1.0_wp + 0.61_wp * q_init(k) )
466	flux = ( ( pt_init(k+1) - pt_init(k-1) ) + &
467	0.61_wp * ( pt_init(k+1) * q_init(k+1) - &
468	pt_init(k-1) * q_init(k-1) ) &
469	) * dd2zu(k)
470	ENDIF
471
472	!
473	!-- Calculate dissipation rate if no prognostic equation is used for
474	!-- dissipation rate
475	IF ( dissipation_1d == 'detering' ) THEN
476	diss1d(k) = c_m*3 e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
477	ELSEIF ( dissipation_1d == 'as_in_3d_model' ) THEN
478	diss1d(k) = ( 0.19_wp + 0.74_wp * l1d_diss(k) / l1d_init(k) &
479	) * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
480	ENDIF
481	!
482	!-- TKE
483	te_e(k) = km1d(k) * ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2&
484	+ ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2&
485	) &
486	- g / pt_0 * kh1d(k) * flux &
487	+ ( &
488	kmzp * ( e1d(k+1) - e1d(k) ) * ddzu(k+1) &
489	- kmzm * ( e1d(k) - e1d(k-1) ) * ddzu(k) &
490	) * ddzw(k) &
491	- diss1d(k)
492
493	IF ( dissipation_1d == 'prognostic' ) THEN
494	!
495	!-- dissipation rate
496	te_diss(k) = km1d(k) * &
497	( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2 &
498	+ ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2 &
499	) * c_1 * c_mu*0.75 / c_h f / us1d &
500	* SQRT(e1d(k)) &
501	- g / pt_0 * kh1d(k) * flux * c_3 &
502	* diss1d(k) / ( e1d(k) + 1.0E-20_wp ) &
503	+ ( kmzp * ( diss1d(k+1) - diss1d(k) ) &
504	* ddzu(k+1) &
505	- kmzm * ( diss1d(k) - diss1d(k-1) ) &
506	* ddzu(k) &
507	) * ddzw(k) / sig_diss &
508	- c_2 * diss1d(k)**2 / ( e1d(k) + 1.0E-20_wp )
509
510	ENDIF
511
512	ENDDO
513	ENDIF
514
515	!
516	!-- Tendency terms at the top of the constant-flux layer.
517	!-- Finite differences of the momentum fluxes are computed using half the
518	!-- normal grid length (2.0*ddzw(k)) for the sake of enhanced accuracy
519	IF ( constant_flux_layer ) THEN
520
521	k = nzb+1
522	kmzm = 0.5_wp * ( km1d(k-1) + km1d(k) )
523	kmzp = 0.5_wp * ( km1d(k) + km1d(k+1) )
524	IF ( .NOT. humidity ) THEN
525	pt_0 = pt_init(k)
526	flux = ( pt_init(k+1)-pt_init(k-1) ) * dd2zu(k)
527	ELSE
528	pt_0 = pt_init(k) * ( 1.0_wp + 0.61_wp * q_init(k) )
529	flux = ( ( pt_init(k+1) - pt_init(k-1) ) + &
530	0.61_wp * ( pt_init(k+1) * q_init(k+1) - &
531	pt_init(k-1) * q_init(k-1) ) &
532	) * dd2zu(k)
533	ENDIF
534
535	!
536	!-- Calculate dissipation rate if no prognostic equation is used for
537	!-- dissipation rate
538	IF ( dissipation_1d == 'detering' ) THEN
539	diss1d(k) = c_m*3 e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
540	ELSEIF ( dissipation_1d == 'as_in_3d_model' ) THEN
541	diss1d(k) = ( 0.19_wp + 0.74_wp * l1d_diss(k) / l1d_init(k) ) &
542	* e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
543	ENDIF
544
545	!
546	!-- u-component
547	te_u(k) = f * ( v1d(k) - vg(k) ) + ( &
548	kmzp * ( u1d(k+1) - u1d(k) ) * ddzu(k+1) + usws1d &
549	) * 2.0_wp * ddzw(k)
550	!
551	!-- v-component
552	te_v(k) = -f * ( u1d(k) - ug(k) ) + ( &
553	kmzp * ( v1d(k+1) - v1d(k) ) * ddzu(k+1) + vsws1d &
554	) * 2.0_wp * ddzw(k)
555	!
556	!-- TKE
557	IF ( .NOT. dissipation_1d == 'prognostic' ) THEN
558	te_e(k) = km1d(k) * ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2&
559	+ ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2&
560	) &
561	- g / pt_0 * kh1d(k) * flux &
562	+ ( &
563	kmzp * ( e1d(k+1) - e1d(k) ) * ddzu(k+1) &
564	- kmzm * ( e1d(k) - e1d(k-1) ) * ddzu(k) &
565	) * ddzw(k) &
566	- diss1d(k)
567	ENDIF
568
569	ENDIF
570
571	!
572	!-- Prognostic equations for all 1D variables
573	DO k = nzb+1, nzt
574
575	u1d_p(k) = u1d(k) + dt_1d * ( tsc(2) * te_u(k) + &
576	tsc(3) * te_um(k) )
577	v1d_p(k) = v1d(k) + dt_1d * ( tsc(2) * te_v(k) + &
578	tsc(3) * te_vm(k) )
579
580	ENDDO
581	IF ( .NOT. constant_diffusion ) THEN
582
583	DO k = nzb+1, nzt
584	e1d_p(k) = e1d(k) + dt_1d * ( tsc(2) * te_e(k) + &
585	tsc(3) * te_em(k) )
586	ENDDO
587
588	IF ( dissipation_1d == 'prognostic' ) THEN
589	DO k = nzb_diff, nzt
590	diss1d_p(k) = diss1d(k) + dt_1d * ( tsc(2) * te_diss(k) + &
591	tsc(3) * te_dissm(k) )
592	ENDDO
593	ENDIF
594	!
595	!-- Eliminate negative TKE values, which can result from the
596	!-- integration due to numerical inaccuracies. In such cases the TKE
597	!-- value is reduced to 10 percent of its old value.
598	WHERE ( e1d_p < 0.0_wp ) e1d_p = 0.1_wp * e1d
599	ENDIF
600
601	!
602	!-- Calculate tendencies for the next Runge-Kutta step
603	IF ( timestep_scheme(1:5) == 'runge' ) THEN
604	IF ( intermediate_timestep_count == 1 ) THEN
605
606	DO k = nzb+1, nzt
607	te_um(k) = te_u(k)
608	te_vm(k) = te_v(k)
609	ENDDO
610
611	IF ( .NOT. constant_diffusion ) THEN
612	DO k = nzb+1, nzt
613	te_em(k) = te_e(k)
614	ENDDO
615	IF ( dissipation_1d == 'prognostic' ) THEN
616	DO k = nzb+1, nzt
617	te_dissm(k) = te_diss(k)
618	ENDDO
619	ENDIF
620	ENDIF
621
622	ELSEIF ( intermediate_timestep_count < &
623	intermediate_timestep_count_max ) THEN
624
625	DO k = nzb+1, nzt
626	te_um(k) = -9.5625_wp * te_u(k) + 5.3125_wp * te_um(k)
627	te_vm(k) = -9.5625_wp * te_v(k) + 5.3125_wp * te_vm(k)
628	ENDDO
629
630	IF ( .NOT. constant_diffusion ) THEN
631	DO k = nzb+1, nzt
632	te_em(k) = -9.5625_wp * te_e(k) + 5.3125_wp * te_em(k)
633	ENDDO
634	IF ( dissipation_1d == 'prognostic' ) THEN
635	DO k = nzb+1, nzt
636	te_dissm(k) = -9.5625_wp * te_diss(k) + 5.3125_wp * te_dissm(k)
637	ENDDO
638	ENDIF
639	ENDIF
640
641	ENDIF
642	ENDIF
643
644
645	!
646	!-- Boundary conditions for the prognostic variables.
647	!-- At the top boundary (nzt+1) u, v, e, and diss keep their initial
648	!-- values (ug(nzt+1), vg(nzt+1), 0, 0).
649	!-- At the bottom boundary, Dirichlet condition is used for u and v (0)
650	!-- and Neumann condition for e and diss (e(nzb)=e(nzb+1)).
651	u1d_p(nzb) = 0.0_wp
652	v1d_p(nzb) = 0.0_wp
653
654	!
655	!-- Swap the time levels in preparation for the next time step.
656	u1d = u1d_p
657	v1d = v1d_p
658	IF ( .NOT. constant_diffusion ) THEN
659	e1d = e1d_p
660	IF ( dissipation_1d == 'prognostic' ) THEN
661	diss1d = diss1d_p
662	ENDIF
663	ENDIF
664
665	!
666	!-- Compute diffusion quantities
667	IF ( .NOT. constant_diffusion ) THEN
668
669	!
670	!-- First compute the vertical fluxes in the constant-flux layer
671	IF ( constant_flux_layer ) THEN
672	!
673	!-- Compute theta* using Rif numbers of the previous time step
674	IF ( rif1d(nzb+1) >= 0.0_wp ) THEN
675	!
676	!-- Stable stratification
677	ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) / &
678	( LOG( zu(nzb+1) / z0h1d ) + 5.0_wp * rif1d(nzb+1) * &
679	( zu(nzb+1) - z0h1d ) / zu(nzb+1) &
680	)
681	ELSE
682	!
683	!-- Unstable stratification
684	a = SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) )
685	b = SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) / &
686	zu(nzb+1) * z0h1d )
687
688	ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) / &
689	LOG( (a-1.0_wp) / (a+1.0_wp) * &
690	(b+1.0_wp) / (b-1.0_wp) )
691	ENDIF
692
693	ENDIF ! constant_flux_layer
694
695	!
696	!-- Compute the Richardson-flux numbers,
697	!-- first at the top of the constant-flux layer using u* of the
698	!-- previous time step (+1E-30, if u* = 0), then in the remaining area.
699	!-- There the rif-numbers of the previous time step are used.
700
701	IF ( constant_flux_layer ) THEN
702	IF ( .NOT. humidity ) THEN
703	pt_0 = pt_init(nzb+1)
704	flux = ts1d
705	ELSE
706	pt_0 = pt_init(nzb+1) * ( 1.0_wp + 0.61_wp * q_init(nzb+1) )
707	flux = ts1d + 0.61_wp * pt_init(k) * qs1d
708	ENDIF
709	rif1d(nzb+1) = zu(nzb+1) * kappa * g * flux / &
710	( pt_0 * ( us1d**2 + 1E-30_wp ) )
711	ENDIF
712
713	DO k = nzb_diff, nzt
714	IF ( .NOT. humidity ) THEN
715	pt_0 = pt_init(k)
716	flux = ( pt_init(k+1) - pt_init(k-1) ) * dd2zu(k)
717	ELSE
718	pt_0 = pt_init(k) * ( 1.0_wp + 0.61_wp * q_init(k) )
719	flux = ( ( pt_init(k+1) - pt_init(k-1) ) &
720	+ 0.61_wp &
721	* ( pt_init(k+1) * q_init(k+1) &
722	- pt_init(k-1) * q_init(k-1) ) &
723	) * dd2zu(k)
724	ENDIF
725	IF ( rif1d(k) >= 0.0_wp ) THEN
726	rif1d(k) = g / pt_0 * flux / &
727	( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2 &
728	+ ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2 &
729	+ 1E-30_wp &
730	)
731	ELSE
732	rif1d(k) = g / pt_0 * flux / &
733	( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2 &
734	+ ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2 &
735	+ 1E-30_wp &
736	) * ( 1.0_wp - 16.0_wp * rif1d(k) )**0.25_wp
737	ENDIF
738	ENDDO
739	!
740	!-- Richardson-numbers must remain restricted to a realistic value
741	!-- range. It is exceeded excessively for very small velocities
742	!-- (u,v --> 0).
743	WHERE ( rif1d < -5.0_wp ) rif1d = -5.0_wp
744	WHERE ( rif1d > 1.0_wp ) rif1d = 1.0_wp
745
746	!
747	!-- Compute u* from the absolute velocity value
748	IF ( constant_flux_layer ) THEN
749	uv_total = SQRT( u1d(nzb+1)2 + v1d(nzb+1)2 )
750
751	IF ( rif1d(nzb+1) >= 0.0_wp ) THEN
752	!
753	!-- Stable stratification
754	us1d = kappa * uv_total / ( &
755	LOG( zu(nzb+1) / z01d ) + 5.0_wp * rif1d(nzb+1) * &
756	( zu(nzb+1) - z01d ) / zu(nzb+1) &
757	)
758	ELSE
759	!
760	!-- Unstable stratification
761	a = 1.0_wp / SQRT( SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) ) )
762	b = 1.0_wp / SQRT( SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) / &
763	zu(nzb+1) * z01d ) )
764	us1d = kappa * uv_total / ( &
765	LOG( (1.0_wp+b) / (1.0_wp-b) * (1.0_wp-a) / &
766	(1.0_wp+a) ) + &
767	2.0_wp * ( ATAN( b ) - ATAN( a ) ) &
768	)
769	ENDIF
770
771	!
772	!-- Compute the momentum fluxes for the diffusion terms
773	usws1d = - u1d(nzb+1) / uv_total * us1d**2
774	vsws1d = - v1d(nzb+1) / uv_total * us1d**2
775
776	!
777	!-- Boundary condition for the turbulent kinetic energy and
778	!-- dissipation rate at the top of the constant-flux layer.
779	!-- Additional Neumann condition de/dz = 0 at nzb is set to ensure
780	!-- compatibility with the 3D model.
781	IF ( ibc_e_b == 2 ) THEN
782	e1d(nzb+1) = ( us1d / c_m )**2
783	ENDIF
784	IF ( dissipation_1d == 'prognostic' ) THEN
785	e1d(nzb+1) = us1d**2 / SQRT( c_mu )
786	diss1d(nzb+1) = us1d*3 / ( kappa zu(nzb+1) )
787	diss1d(nzb) = diss1d(nzb+1)
788	ENDIF
789	e1d(nzb) = e1d(nzb+1)
790
791	IF ( humidity ) THEN
792	!
793	!-- Compute q*
794	IF ( rif1d(nzb+1) >= 0.0_wp ) THEN
795	!
796	!-- Stable stratification
797	qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) / &
798	( LOG( zu(nzb+1) / z0h1d ) + 5.0_wp * rif1d(nzb+1) * &
799	( zu(nzb+1) - z0h1d ) / zu(nzb+1) &
800	)
801	ELSE
802	!
803	!-- Unstable stratification
804	a = SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) )
805	b = SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) / &
806	zu(nzb+1) * z0h1d )
807	qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) / &
808	LOG( (a-1.0_wp) / (a+1.0_wp) * &
809	(b+1.0_wp) / (b-1.0_wp) )
810	ENDIF
811	ELSE
812	qs1d = 0.0_wp
813	ENDIF
814
815	ENDIF ! constant_flux_layer
816
817	!
818	!-- Compute the diabatic mixing length. The unstable stratification
819	!-- must not be considered for l1d (km1d) as it is already considered
820	!-- in the dissipation of TKE via l1d_diss. Otherwise, km1d would be
821	!-- too large.
822	IF ( mixing_length_1d == 'blackadar' ) THEN
823	DO k = nzb+1, nzt
824	IF ( rif1d(k) >= 0.0_wp ) THEN
825	l1d(k) = l1d_init(k) / ( 1.0_wp + 5.0_wp * rif1d(k) )
826	l1d_diss(k) = l1d(k)
827	ELSE
828	l1d(k) = l1d_init(k)
829	l1d_diss(k) = l1d_init(k) * &
830	SQRT( 1.0_wp - 16.0_wp * rif1d(k) )
831	ENDIF
832	ENDDO
833	ELSEIF ( mixing_length_1d == 'as_in_3d_model' ) THEN
834	DO k = nzb+1, nzt
835	dpt_dz = ( pt_init(k+1) - pt_init(k-1) ) * dd2zu(k)
836	IF ( dpt_dz > 0.0_wp ) THEN
837	l_stable = 0.76_wp * SQRT( e1d(k) ) / &
838	SQRT( g / pt_init(k) * dpt_dz ) + 1E-5_wp
839	ELSE
840	l_stable = l1d_init(k)
841	ENDIF
842	l1d(k) = MIN( l1d_init(k), l_stable )
843	l1d_diss(k) = l1d(k)
844	ENDDO
845	ENDIF
846
847	!
848	!-- Compute the diffusion coefficients for momentum via the
849	!-- corresponding Prandtl-layer relationship and according to
850	!-- Prandtl-Kolmogorov, respectively
851	IF ( constant_flux_layer ) THEN
852	IF ( rif1d(nzb+1) >= 0.0_wp ) THEN
853	km1d(nzb+1) = us1d * kappa * zu(nzb+1) / &
854	( 1.0_wp + 5.0_wp * rif1d(nzb+1) )
855	ELSE
856	km1d(nzb+1) = us1d * kappa * zu(nzb+1) * &
857	( 1.0_wp - 16.0_wp * rif1d(nzb+1) )**0.25_wp
858	ENDIF
859	ENDIF
860	IF ( dissipation_1d == 'prognostic' ) THEN
861	DO k = nzb_diff, nzt
862	km1d(k) = c_mu * e1d(k)**2 / ( diss1d(k) + 1.0E-30_wp )
863	ENDDO
864	ELSE
865	DO k = nzb_diff, nzt
866	km1d(k) = c_m * SQRT( e1d(k) ) * l1d(k)
867	ENDDO
868	ENDIF
869
870	!
871	!-- Add damping layer
872	DO k = damp_level_ind_1d+1, nzt+1
873	km1d(k) = 1.1_wp * km1d(k-1)
874	km1d(k) = MIN( km1d(k), 10.0_wp )
875	ENDDO
876
877	!
878	!-- Compute the diffusion coefficient for heat via the relationship
879	!-- kh = phim / phih * km
880	DO k = nzb+1, nzt
881	IF ( rif1d(k) >= 0.0_wp ) THEN
882	kh1d(k) = km1d(k)
883	ELSE
884	kh1d(k) = km1d(k) * ( 1.0_wp - 16.0_wp * rif1d(k) )**0.25_wp
885	ENDIF
886	ENDDO
887
888	ENDIF ! .NOT. constant_diffusion
889
890	ENDDO ! intermediate step loop
891
892	!
893	!-- Increment simulated time and output times
894	current_timestep_number_1d = current_timestep_number_1d + 1
895	simulated_time_1d = simulated_time_1d + dt_1d
896	simulated_time_chr = time_to_string( simulated_time_1d )
897	time_pr_1d = time_pr_1d + dt_1d
898	time_run_control_1d = time_run_control_1d + dt_1d
899
900	!
901	!-- Determine and print out quantities for run control
902	IF ( time_run_control_1d >= dt_run_control_1d ) THEN
903	CALL run_control_1d
904	time_run_control_1d = time_run_control_1d - dt_run_control_1d
905	ENDIF
906
907	!
908	!-- Profile output on file
909	IF ( time_pr_1d >= dt_pr_1d ) THEN
910	CALL print_1d_model
911	time_pr_1d = time_pr_1d - dt_pr_1d
912	ENDIF
913
914	!
915	!-- Determine size of next time step
916	CALL timestep_1d
917
918	ENDDO ! time loop
919
920
921	END SUBROUTINE time_integration_1d
922
923
924	!------------------------------------------------------------------------------!
925	! Description:
926	! ------------
927	!> Compute and print out quantities for run control of the 1D model.
928	!------------------------------------------------------------------------------!
929
930	SUBROUTINE run_control_1d
931
932
933	USE constants, &
934	ONLY: pi
935
936	IMPLICIT NONE
937
938	INTEGER(iwp) :: k !< loop index
939
940	REAL(wp) :: alpha !< angle of wind vector at top of constant-flux layer
941	REAL(wp) :: energy !< kinetic energy
942	REAL(wp) :: umax !< maximum of u
943	REAL(wp) :: uv_total !< horizontal wind speed
944	REAL(wp) :: vmax !< maximum of v
945
946	!
947	!-- Output
948	IF ( myid == 0 ) THEN
949	!
950	!-- If necessary, write header
951	IF ( .NOT. run_control_header_1d ) THEN
952	CALL check_open( 15 )
953	WRITE ( 15, 100 )
954	run_control_header_1d = .TRUE.
955	ENDIF
956
957	!
958	!-- Compute control quantities
959	!-- grid level nzp is excluded due to mirror boundary condition
960	umax = 0.0_wp; vmax = 0.0_wp; energy = 0.0_wp
961	DO k = nzb+1, nzt+1
962	umax = MAX( ABS( umax ), ABS( u1d(k) ) )
963	vmax = MAX( ABS( vmax ), ABS( v1d(k) ) )
964	energy = energy + 0.5_wp * ( u1d(k)2 + v1d(k)2 )
965	ENDDO
966	energy = energy / REAL( nzt - nzb + 1, KIND=wp )
967
968	uv_total = SQRT( u1d(nzb+1)2 + v1d(nzb+1)2 )
969	IF ( ABS( v1d(nzb+1) ) < 1.0E-5_wp ) THEN
970	alpha = ACOS( SIGN( 1.0_wp , u1d(nzb+1) ) )
971	ELSE
972	alpha = ACOS( u1d(nzb+1) / uv_total )
973	IF ( v1d(nzb+1) <= 0.0_wp ) alpha = 2.0_wp * pi - alpha
974	ENDIF
975	alpha = alpha / ( 2.0_wp * pi ) * 360.0_wp
976
977	WRITE ( 15, 101 ) current_timestep_number_1d, simulated_time_chr, &
978	dt_1d, umax, vmax, us1d, alpha, energy
979	!
980	!-- Write buffer contents to disc immediately
981	FLUSH( 15 )
982
983	ENDIF
984
985	!
986	!-- formats
987	100 FORMAT (///'1D run control output:'/ &
988	&'------------------------------'// &
989	&'ITER. HH:MM:SS DT UMAX VMAX U* ALPHA ENERG.'/ &
990	&'-------------------------------------------------------------')
991	101 FORMAT (I5,2X,A9,1X,F6.2,2X,F6.2,1X,F6.2,1X,F6.3,2X,F5.1,2X,F7.2)
992
993
994	END SUBROUTINE run_control_1d
995
996
997
998	!------------------------------------------------------------------------------!
999	! Description:
1000	! ------------
1001	!> Compute the time step w.r.t. the diffusion criterion
1002	!------------------------------------------------------------------------------!
1003
1004	SUBROUTINE timestep_1d
1005
1006	IMPLICIT NONE
1007
1008	INTEGER(iwp) :: k !< loop index
1009
1010	REAL(wp) :: dt_diff !< time step accorind to diffusion criterion
1011	REAL(wp) :: fac !< factor of criterion
1012	REAL(wp) :: value !< auxiliary variable
1013
1014	!
1015	!-- Compute the currently feasible time step according to the diffusion
1016	!-- criterion. At nzb+1 the half grid length is used.
1017	fac = 0.125 !0.35_wp !### changed from 0.35
1018	dt_diff = dt_max_1d
1019	DO k = nzb+2, nzt
1020	value = fac * dzu(k) * dzu(k) / ( km1d(k) + 1E-20_wp )
1021	dt_diff = MIN( value, dt_diff )
1022	ENDDO
1023	value = fac * zu(nzb+1) * zu(nzb+1) / ( km1d(nzb+1) + 1E-20_wp )
1024	dt_1d = MIN( value, dt_diff )
1025
1026	!
1027	!-- Set flag when the time step becomes too small
1028	IF ( dt_1d < ( 0.00001_wp * dt_max_1d ) ) THEN
1029	stop_dt_1d = .TRUE.
1030
1031	WRITE( message_string, * ) 'timestep has exceeded the lower limit &', &
1032	'dt_1d = ',dt_1d,' s simulation stopped!'
1033	CALL message( 'timestep_1d', 'PA0192', 1, 2, 0, 6, 0 )
1034
1035	ENDIF
1036
1037	END SUBROUTINE timestep_1d
1038
1039
1040
1041	!------------------------------------------------------------------------------!
1042	! Description:
1043	! ------------
1044	!> List output of profiles from the 1D-model
1045	!------------------------------------------------------------------------------!
1046
1047	SUBROUTINE print_1d_model
1048
1049	IMPLICIT NONE
1050
1051	INTEGER(iwp) :: k !< loop parameter
1052
1053	LOGICAL, SAVE :: write_first = .TRUE. !< flag for writing header
1054
1055
1056	IF ( myid == 0 ) THEN
1057	!
1058	!-- Open list output file for profiles from the 1D-model
1059	CALL check_open( 17 )
1060
1061	!
1062	!-- Write Header
1063	IF ( write_first ) THEN
1064	WRITE ( 17, 100 ) TRIM( run_description_header )
1065	write_first = .FALSE.
1066	ENDIF
1067
1068	!
1069	!-- Write the values
1070	WRITE ( 17, 104 ) TRIM( simulated_time_chr )
1071	WRITE ( 17, 101 )
1072	WRITE ( 17, 102 )
1073	WRITE ( 17, 101 )
1074	DO k = nzt+1, nzb, -1
1075	WRITE ( 17, 103) k, zu(k), u1d(k), v1d(k), pt_init(k), e1d(k), &
1076	rif1d(k), km1d(k), kh1d(k), l1d(k), diss1d(k)
1077	ENDDO
1078	WRITE ( 17, 101 )
1079	WRITE ( 17, 102 )
1080	WRITE ( 17, 101 )
1081
1082	!
1083	!-- Write buffer contents to disc immediately
1084	FLUSH( 17 )
1085
1086	ENDIF
1087
1088	!
1089	!-- Formats
1090	100 FORMAT ('# ',A/'#',10('-')/'# 1d-model profiles')
1091	104 FORMAT (//'# Time: ',A)
1092	101 FORMAT ('#',111('-'))
1093	102 FORMAT ('# k zu u v pt e ', &
1094	'rif Km Kh l diss ')
1095	103 FORMAT (1X,I4,1X,F7.1,9(1X,E10.3))
1096
1097
1098	END SUBROUTINE print_1d_model
1099
1100
1101	END MODULE

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |