Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |