Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |