Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |