Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |