Home

Context Navigation

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

Last change on this file since 2372 was 2339, checked in by gronemeier, 7 years ago
corrected timestamp in header
Property svn:keywords set to `Id`
File size: 37.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |