Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |