Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |