Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |