Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |