Home

Context Navigation

source: palm/trunk/SOURCE/model_1d_mod.f90 @ 4181

Last change on this file since 4181 was 4180, checked in by scharf, 5 years ago
removed comments in 'Former revisions' section that are older than 01.01.2019
Property svn:keywords set to `Id`
File size: 40.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |