Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |