Home

Context Navigation

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

Last change on this file since 54 was 46, checked in by raasch, 18 years ago
updating parts of Marcus changes
Property svn:keywords set to `Id`
File size: 33.7 KB

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