Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |