Home

Context Navigation

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

Last change on this file since 1683 was 1683, checked in by knoop, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 35.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |