Home

Context Navigation

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

Last change on this file since 16 was 4, checked in by raasch, 18 years ago
Id keyword set as property for all *.f90 files
Property svn:keywords set to `Id`
File size: 33.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |