Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |