Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |