Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |