Home

Context Navigation

source: palm/trunk/SOURCE/init_1d_model.f90 @ 1806

Last change on this file since 1806 was 1710, checked in by maronga, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 36.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |