Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |