Home

Context Navigation

model_1d_mod.f90 @ 4894

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

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

Context Navigation

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

Download in other formats: