Home

Context Navigation

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

Last change on this file since 1576 was 1354, checked in by heinze, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 35.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |