Home

Context Navigation

init_1d_model.f90 @ 419

Last change on this file since 419 was 392, checked in by raasch, 15 years ago

New:
---

Adapted for machine lck
(mrun, mbuild, subjob)

bc_lr/bc_ns in most subroutines replaced by LOGICAL variables bc_lr_cyc,
bc_ns_cyc for speed optimization
(check_parameters, diffusion_u, diffusion_v, diffusion_w, modules)

Additional timestep criterion in case of simulations with plant canopy (timestep)

Check for illegal entries in section_xy|xz|yz that exceed nz+1|ny+1|nx+1
(check_parameters)

Clipping of dvrp output implemented. Default colourtable for particles
implemented, particle attributes (color, dvrp_size) can be set with new
parameters particle_color, particle_dvrpsize, color_interval,
dvrpsize_interval (init_dvrp, data_output_dvrp, modules, user_data_output_dvrp).
Slicer attributes (dvrp) are set with new routine set_slicer_attributes_dvrp
and are controlled with existing parameters slicer_range_limits.
(set_slicer_attributes_dvrp)

Ocean atmosphere coupling allows to use independent precursor runs in order
to account for different spin-up times. The time when coupling has to be
started is given by new inipar parameter coupling_start_time. The precursor
ocean run has to be started using new mrun option "-y" in order to add
appendix "_O" to all output files.
(check_for_restart, check_parameters, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, header, init_coupling, modules, mrun,
parin, read_var_list, surface_coupler, time_integration, write_var_list)

Polygon reduction for topography and ground plate isosurface. Reduction level
for buildings can be chosen with parameter cluster_size. (init_dvrp)

External pressure gradient (check_parameters, header, init_3d_model, modules,
parin, prognostic_equations, read_var_list, write_var_list)

New topography case 'single_street_canyon' (header, init_grid, modules, parin,
read_var_list, user_check_parameters, user_header, user_init_grid, write_var_list)

Option to predefine a target bulk velocity for conserve_volume_flow
(check_parameters, header, init_3d_model, modules, parin, read_var_list,
write_var_list)

Option for user defined 2D data output in xy cross sections at z=nzb+1
(data_output_2d, user_data_output_2d)

xy cross section output of surface heatfluxes (latent, sensible)
(average_3d_data, check_parameters, data_output_2d, modules, read_3d_binary,
sum_up_3d_data, write_3d_binary)

average_3d_data, check_for_restart, check_parameters, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, init_coupling, init_dvrp, init_grid, init_3d_model, header, mbuild, modules, mrun, package_parin, parin, prognostic_equations, read_3d_binary, read_var_list, subjob, surface_coupler, timestep, time_integration, user_check_parameters, user_data_output_2d, user_data_output_dvrp, user_header, user_init_grid, write_3d_binary, write_var_list

New: set_particle_attributes, set_slicer_attributes_dvrp

Changed:

lcmuk changed to lc to avoid problems with Intel compiler on sgi-ice
(poisfft)

For extended NetCDF files, the updated title attribute includes an update of
time_average_text where appropriate. (netcdf)

In case of restart runs without extension, initial profiles are not written
to NetCDF-file anymore. (data_output_profiles, modules, read_var_list, write_var_list)

Small change in formatting of the message handling routine concering the output in the
job protocoll. (message)

initializing_actions='read_data_for_recycling' renamed to 'cyclic_fill', now
independent of turbulent_inflow (check_parameters, header, init_3d_model)

2 NetCDF error numbers changed. (data_output_3d)

A Link to the website appendix_a.html is printed for further information
about the possible errors. (message)

Temperature gradient criterion for estimating the boundary layer height
replaced by the gradient criterion of Sullivan et al. (1998). (flow_statistics)

NetCDF unit attribute in timeseries output in case of statistic regions added
(netcdf)

Output of NetCDF messages with aid of message handling routine.
(check_open, close_file, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, netcdf, output_particles_netcdf)

Output of messages replaced by message handling routine.
(advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart,
check_open, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp,
data_output_profiles, data_output_spectra, fft_xy, flow_statistics, header,
init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid,
netcdf, parin, plant_canopy_model, poisfft_hybrid, poismg, read_3d_binary,
read_var_list, surface_coupler, temperton_fft, timestep, user_actions,
user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy,
user_parin, user_read_restart_data, user_spectra )

Maximum number of tails is calculated from maximum number of particles and
skip_particles_for_tail (init_particles)

Value of vertical_particle_advection may differ for each particle group
(advec_particles, header, modules)

First constant in array den also defined as type double. (eqn_state_seawater)

Parameter dvrp_psize moved from particles_par to dvrp_graphics_par. (package_parin)

topography_grid_convention moved from userpar to inipar (check_parameters,
header, parin, read_var_list, user_check_parameters, user_header,
user_init_grid, user_parin, write_var_list)

Default value of grid_matching changed to strict.

Adjustments for runs on lcxt4 (necessary due to an software update on CRAY) and
for coupled runs on ibmy (mrun, subjob)

advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart, check_open, check_parameters, close_file, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, eqn_state_seawater, fft_xy, flow_statistics, header, init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid, message, mrun, netcdf, output_particles_netcdf, package_parin, parin, plant_canopy_model, poisfft, poisfft_hybrid, poismg, read_3d_binary, read_var_list, sort_particles, subjob, user_check_parameters, user_header, user_init_grid, user_parin, surface_coupler, temperton_fft, timestep, user_actions, user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy, user_parin, user_read_restart_data, user_spectra, write_var_list

Errors:

Bugfix: Initial hydrostatic pressure profile in case of ocean runs is now
calculated in 5 iteration steps. (init_ocean)

Bugfix: wrong sign in buoyancy production of ocean part in case of not using
the reference density (only in 3D routine production_e) (production_e)

Bugfix: output of averaged 2d/3d quantities requires that an avaraging
interval has been set, respective error message is included (check_parameters)

Bugfix: Output on unit 14 only if requested by write_binary.
(user_last_actions)

Bugfix to avoid zero division by km_neutral (production_e)

Bugfix for extended NetCDF files: In order to avoid 'data mode' errors if
updated attributes are larger than their original size, NF90_PUT_ATT is called
in 'define mode' enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a
possible performance loss; an alternative strategy would be to ensure equal
attribute size in a job chain. (netcdf)

Bugfix: correction of initial volume flow for non-flat topography (init_3d_model)
Bugfix: zero initialization of arrays within buildings for 'cyclic_fill' (init_3d_model)

Bugfix: to_be_resorted => s_av for time-averaged scalars (data_output_2d, data_output_3d)

Bugfix: error in formatting the output (message)

Bugfix: avoid that ngp_2dh_s_inner becomes zero (init_3_model)

Typographical error: unit of wpt in dots_unit (modules)

Bugfix: error in check, if particles moved further than one subdomain length.
This check must not be applied for newly released particles. (advec_particles)

Bugfix: several tail counters are initialized, particle_tail_coordinates is
only written to file if its third index is > 0, arrays for tails are allocated
with a minimum size of 10 tails if there is no tail initially (init_particles,
advec_particles)

Bugfix: pressure included for profile output (check_parameters)

Bugfix: Type of count and count_rate changed to default INTEGER on NEC machines
(cpu_log)

Bugfix: output if particle time series only if particle advection is switched
on. (time_integration)

Bugfix: qsws was calculated in case of constant heatflux = .FALSE. (prandtl_fluxes)

Bugfix: averaging along z is not allowed for 2d quantities (e.g. u* and z0) (data_output_2d)

Typographical errors (netcdf)

If the inversion height calculated by the prerun is zero, inflow_damping_height
must be explicitly specified (init_3d_model)

Small bugfix concerning 3d 64bit netcdf output format (header)

Bugfix: dt_fixed removed from the restart file, because otherwise, no change
from a fixed to a variable timestep would be possible in restart runs.
(read_var_list, write_var_list)

Bugfix: initial setting of time_coupling in coupled restart runs (time_integration)

advec_particles, check_parameters, cpu_log, data_output_2d, data_output_3d, header, init_3d_model, init_particles, init_ocean, modules, netcdf, prandtl_fluxes, production_e, read_var_list, time_integration, user_last_actions, write_var_list

Property svn:keywords set to Id

File size: 34.1 KB

Line
1	SUBROUTINE init_1d_model
2
3	!------------------------------------------------------------------------------!
4	! Current revisions:
5	! -----------------
6	!
7	!
8	! Former revisions:
9	! -----------------
10	! $Id: init_1d_model.f90 392 2009-09-24 10:39:14Z heinze $
11	!
12	! 254 2009-03-05 15:33:42Z heinze
13	! Output of messages replaced by message handling routine.
14	!
15	! 184 2008-08-04 15:53:39Z letzel
16	! provisional solution for run_control_1d output: add 'CALL check_open( 15 )'
17	!
18	! 135 2007-11-22 12:24:23Z raasch
19	! Bugfix: absolute value of f must be used when calculating the Blackadar
20	! mixing length
21	!
22	! 82 2007-04-16 15:40:52Z raasch
23	! Preprocessor strings for different linux clusters changed to "lc",
24	! routine local_flush is used for buffer flushing
25	!
26	! 75 2007-03-22 09:54:05Z raasch
27	! Bugfix: preset of tendencies te_em, te_um, te_vm,
28	! moisture renamed humidity
29	!
30	! RCS Log replace by Id keyword, revision history cleaned up
31	!
32	! Revision 1.21 2006/06/02 15:19:57 raasch
33	! cpp-directives extended for lctit
34	!
35	! Revision 1.1 1998/03/09 16:22:10 raasch
36	! Initial revision
37	!
38	!
39	! Description:
40	! ------------
41	! 1D-model to initialize the 3D-arrays.
42	! The temperature profile is set as steady and a corresponding steady solution
43	! of the wind profile is being computed.
44	! All subroutines required can be found within this file.
45	!------------------------------------------------------------------------------!
46
47	USE arrays_3d
48	USE indices
49	USE model_1d
50	USE control_parameters
51
52	IMPLICIT NONE
53
54	CHARACTER (LEN=9) :: time_to_string
55	INTEGER :: k
56	REAL :: lambda
57
58	!
59	!-- Allocate required 1D-arrays
60	ALLOCATE( e1d(nzb:nzt+1), e1d_m(nzb:nzt+1), e1d_p(nzb:nzt+1), &
61	kh1d(nzb:nzt+1), kh1d_m(nzb:nzt+1), km1d(nzb:nzt+1), &
62	km1d_m(nzb:nzt+1), l_black(nzb:nzt+1), l1d(nzb:nzt+1), &
63	l1d_m(nzb:nzt+1), rif1d(nzb:nzt+1), te_e(nzb:nzt+1), &
64	te_em(nzb:nzt+1), te_u(nzb:nzt+1), te_um(nzb:nzt+1), &
65	te_v(nzb:nzt+1), te_vm(nzb:nzt+1), u1d(nzb:nzt+1), &
66	u1d_m(nzb:nzt+1), u1d_p(nzb:nzt+1), v1d(nzb:nzt+1), &
67	v1d_m(nzb:nzt+1), v1d_p(nzb:nzt+1) )
68
69	!
70	!-- Initialize arrays
71	IF ( constant_diffusion ) THEN
72	km1d = km_constant
73	km1d_m = km_constant
74	kh1d = km_constant / prandtl_number
75	kh1d_m = km_constant / prandtl_number
76	ELSE
77	e1d = 0.0; e1d_m = 0.0; e1d_p = 0.0
78	kh1d = 0.0; kh1d_m = 0.0; km1d = 0.0; km1d_m = 0.0
79	rif1d = 0.0
80	!
81	!-- Compute the mixing length
82	l_black(nzb) = 0.0
83
84	IF ( TRIM( mixing_length_1d ) == 'blackadar' ) THEN
85	!
86	!-- Blackadar mixing length
87	IF ( f /= 0.0 ) THEN
88	lambda = 2.7E-4 * SQRT( ug(nzt+1)2 + vg(nzt+1)2 ) / &
89	ABS( f ) + 1E-10
90	ELSE
91	lambda = 30.0
92	ENDIF
93
94	DO k = nzb+1, nzt+1
95	l_black(k) = kappa * zu(k) / ( 1.0 + kappa * zu(k) / lambda )
96	ENDDO
97
98	ELSEIF ( TRIM( mixing_length_1d ) == 'as_in_3d_model' ) THEN
99	!
100	!-- Use the same mixing length as in 3D model
101	l_black(1:nzt) = l_grid
102	l_black(nzt+1) = l_black(nzt)
103
104	ENDIF
105
106	!
107	!-- Adjust mixing length to the prandtl mixing length (within the prandtl
108	!-- layer)
109	IF ( adjust_mixing_length .AND. prandtl_layer ) THEN
110	k = nzb+1
111	l_black(k) = MIN( l_black(k), kappa * zu(k) )
112	ENDIF
113	ENDIF
114	l1d = l_black
115	l1d_m = l_black
116	u1d = u_init
117	u1d_m = u_init
118	u1d_p = u_init
119	v1d = v_init
120	v1d_m = v_init
121	v1d_p = v_init
122
123	!
124	!-- Set initial horizontal velocities at the lowest grid levels to a very small
125	!-- value in order to avoid too small time steps caused by the diffusion limit
126	!-- in the initial phase of a run (at k=1, dz/2 occurs in the limiting formula!)
127	u1d(0:1) = 0.1
128	u1d_m(0:1) = 0.1
129	u1d_p(0:1) = 0.1
130	v1d(0:1) = 0.1
131	v1d_m(0:1) = 0.1
132	v1d_p(0:1) = 0.1
133
134	!
135	!-- For u, theta and the momentum fluxes plausible values are set
136	IF ( prandtl_layer ) THEN
137	us1d = 0.1 ! without initial friction the flow would not change
138	ELSE
139	e1d(nzb+1) = 1.0
140	km1d(nzb+1) = 1.0
141	us1d = 0.0
142	ENDIF
143	ts1d = 0.0
144	usws1d = 0.0; usws1d_m = 0.0
145	vsws1d = 0.0; vsws1d_m = 0.0
146	z01d = roughness_length
147	IF ( humidity .OR. passive_scalar ) qs1d = 0.0
148
149	!
150	!-- Tendencies must be preset in order to avoid runtime errors within the
151	!-- first Runge-Kutta step
152	te_em = 0.0
153	te_um = 0.0
154	te_vm = 0.0
155
156	!
157	!-- Set start time in hh:mm:ss - format
158	simulated_time_chr = time_to_string( simulated_time_1d )
159
160	!
161	!-- Integrate the 1D-model equations using the leap-frog scheme
162	CALL time_integration_1d
163
164
165	END SUBROUTINE init_1d_model
166
167
168
169	SUBROUTINE time_integration_1d
170
171	!------------------------------------------------------------------------------!
172	! Description:
173	! ------------
174	! Leap-frog time differencing scheme for the 1D-model.
175	!------------------------------------------------------------------------------!
176
177	USE arrays_3d
178	USE control_parameters
179	USE indices
180	USE model_1d
181	USE pegrid
182
183	IMPLICIT NONE
184
185	CHARACTER (LEN=9) :: time_to_string
186	INTEGER :: k
187	REAL :: a, b, dissipation, dpt_dz, flux, kmzm, kmzp, l_stable, pt_0, &
188	uv_total
189
190	!
191	!-- Determine the time step at the start of a 1D-simulation and
192	!-- determine and printout quantities used for run control
193	CALL timestep_1d
194	CALL run_control_1d
195
196	!
197	!-- Start of time loop
198	DO WHILE ( simulated_time_1d < end_time_1d .AND. .NOT. stop_dt_1d )
199
200	!
201	!-- Depending on the timestep scheme, carry out one or more intermediate
202	!-- timesteps
203
204	intermediate_timestep_count = 0
205	DO WHILE ( intermediate_timestep_count < &
206	intermediate_timestep_count_max )
207
208	intermediate_timestep_count = intermediate_timestep_count + 1
209
210	CALL timestep_scheme_steering
211
212	!
213	!-- Compute all tendency terms. If a Prandtl-layer is simulated, k starts
214	!-- at nzb+2.
215	DO k = nzb_diff, nzt
216
217	kmzm = 0.5 * ( km1d_m(k-1) + km1d_m(k) )
218	kmzp = 0.5 * ( km1d_m(k) + km1d_m(k+1) )
219	!
220	!-- u-component
221	te_u(k) = f * ( v1d(k) - vg(k) ) + ( &
222	kmzp * ( u1d_m(k+1) - u1d_m(k) ) * ddzu(k+1) &
223	- kmzm * ( u1d_m(k) - u1d_m(k-1) ) * ddzu(k) &
224	) * ddzw(k)
225	!
226	!-- v-component
227	te_v(k) = -f * ( u1d(k) - ug(k) ) + ( &
228	kmzp * ( v1d_m(k+1) - v1d_m(k) ) * ddzu(k+1) &
229	- kmzm * ( v1d_m(k) - v1d_m(k-1) ) * ddzu(k) &
230	) * ddzw(k)
231	ENDDO
232	IF ( .NOT. constant_diffusion ) THEN
233	DO k = nzb_diff, nzt
234	!
235	!-- TKE
236	kmzm = 0.5 * ( km1d_m(k-1) + km1d_m(k) )
237	kmzp = 0.5 * ( km1d_m(k) + km1d_m(k+1) )
238	IF ( .NOT. humidity ) THEN
239	pt_0 = pt_init(k)
240	flux = ( pt_init(k+1)-pt_init(k-1) ) * dd2zu(k)
241	ELSE
242	pt_0 = pt_init(k) * ( 1.0 + 0.61 * q_init(k) )
243	flux = ( ( pt_init(k+1) - pt_init(k-1) ) + &
244	0.61 * pt_init(k) * ( q_init(k+1) - q_init(k-1) ) &
245	) * dd2zu(k)
246	ENDIF
247
248	IF ( dissipation_1d == 'detering' ) THEN
249	!
250	!-- According to Detering, c_e=0.064
251	dissipation = 0.064 * e1d_m(k) * SQRT( e1d_m(k) ) / l1d_m(k)
252	ELSEIF ( dissipation_1d == 'as_in_3d_model' ) THEN
253	dissipation = ( 0.19 + 0.74 * l1d_m(k) / l_grid(k) ) &
254	* e1d_m(k) * SQRT( e1d_m(k) ) / l1d_m(k)
255	ENDIF
256
257	te_e(k) = km1d(k) * ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2&
258	+ ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2&
259	) &
260	- g / pt_0 * kh1d(k) * flux &
261	+ ( &
262	kmzp * ( e1d_m(k+1) - e1d_m(k) ) * ddzu(k+1) &
263	- kmzm * ( e1d_m(k) - e1d_m(k-1) ) * ddzu(k) &
264	) * ddzw(k) &
265	- dissipation
266	ENDDO
267	ENDIF
268
269	!
270	!-- Tendency terms at the top of the Prandtl-layer.
271	!-- Finite differences of the momentum fluxes are computed using half the
272	!-- normal grid length (2.0*ddzw(k)) for the sake of enhanced accuracy
273	IF ( prandtl_layer ) THEN
274
275	k = nzb+1
276	kmzm = 0.5 * ( km1d_m(k-1) + km1d_m(k) )
277	kmzp = 0.5 * ( km1d_m(k) + km1d_m(k+1) )
278	IF ( .NOT. humidity ) THEN
279	pt_0 = pt_init(k)
280	flux = ( pt_init(k+1)-pt_init(k-1) ) * dd2zu(k)
281	ELSE
282	pt_0 = pt_init(k) * ( 1.0 + 0.61 * q_init(k) )
283	flux = ( ( pt_init(k+1) - pt_init(k-1) ) + &
284	0.61 * pt_init(k) * ( q_init(k+1) - q_init(k-1) ) &
285	) * dd2zu(k)
286	ENDIF
287
288	IF ( dissipation_1d == 'detering' ) THEN
289	!
290	!-- According to Detering, c_e=0.064
291	dissipation = 0.064 * e1d_m(k) * SQRT( e1d_m(k) ) / l1d_m(k)
292	ELSEIF ( dissipation_1d == 'as_in_3d_model' ) THEN
293	dissipation = ( 0.19 + 0.74 * l1d_m(k) / l_grid(k) ) &
294	* e1d_m(k) * SQRT( e1d_m(k) ) / l1d_m(k)
295	ENDIF
296
297	!
298	!-- u-component
299	te_u(k) = f * ( v1d(k) - vg(k) ) + ( &
300	kmzp * ( u1d_m(k+1) - u1d_m(k) ) * ddzu(k+1) + usws1d_m &
301	) * 2.0 * ddzw(k)
302	!
303	!-- v-component
304	te_v(k) = -f * ( u1d(k) - ug(k) ) + ( &
305	kmzp * ( v1d_m(k+1) - v1d_m(k) ) * ddzu(k+1) + vsws1d_m &
306	) * 2.0 * ddzw(k)
307	!
308	!-- TKE
309	te_e(k) = km1d(k) * ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2 &
310	+ ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2 &
311	) &
312	- g / pt_0 * kh1d(k) * flux &
313	+ ( &
314	kmzp * ( e1d_m(k+1) - e1d_m(k) ) * ddzu(k+1) &
315	- kmzm * ( e1d_m(k) - e1d_m(k-1) ) * ddzu(k) &
316	) * ddzw(k) &
317	- dissipation
318	ENDIF
319
320	!
321	!-- Prognostic equations for all 1D variables
322	DO k = nzb+1, nzt
323
324	u1d_p(k) = ( 1. - tsc(1) ) * u1d_m(k) + &
325	tsc(1) * u1d(k) + dt_1d * ( tsc(2) * te_u(k) + &
326	tsc(3) * te_um(k) )
327	v1d_p(k) = ( 1. - tsc(1) ) * v1d_m(k) + &
328	tsc(1) * v1d(k) + dt_1d * ( tsc(2) * te_v(k) + &
329	tsc(3) * te_vm(k) )
330
331	ENDDO
332	IF ( .NOT. constant_diffusion ) THEN
333	DO k = nzb+1, nzt
334
335	e1d_p(k) = ( 1. - tsc(1) ) * e1d_m(k) + &
336	tsc(1) * e1d(k) + dt_1d * ( tsc(2) * te_e(k) + &
337	tsc(3) * te_em(k) )
338
339	ENDDO
340	!
341	!-- Eliminate negative TKE values, which can result from the
342	!-- integration due to numerical inaccuracies. In such cases the TKE
343	!-- value is reduced to 10 percent of its old value.
344	WHERE ( e1d_p < 0.0 ) e1d_p = 0.1 * e1d
345	ENDIF
346
347	!
348	!-- Calculate tendencies for the next Runge-Kutta step
349	IF ( timestep_scheme(1:5) == 'runge' ) THEN
350	IF ( intermediate_timestep_count == 1 ) THEN
351
352	DO k = nzb+1, nzt
353	te_um(k) = te_u(k)
354	te_vm(k) = te_v(k)
355	ENDDO
356
357	IF ( .NOT. constant_diffusion ) THEN
358	DO k = nzb+1, nzt
359	te_em(k) = te_e(k)
360	ENDDO
361	ENDIF
362
363	ELSEIF ( intermediate_timestep_count < &
364	intermediate_timestep_count_max ) THEN
365
366	DO k = nzb+1, nzt
367	te_um(k) = -9.5625 * te_u(k) + 5.3125 * te_um(k)
368	te_vm(k) = -9.5625 * te_v(k) + 5.3125 * te_vm(k)
369	ENDDO
370
371	IF ( .NOT. constant_diffusion ) THEN
372	DO k = nzb+1, nzt
373	te_em(k) = -9.5625 * te_e(k) + 5.3125 * te_em(k)
374	ENDDO
375	ENDIF
376
377	ENDIF
378	ENDIF
379
380
381	!
382	!-- Boundary conditions for the prognostic variables.
383	!-- At the top boundary (nzt+1) u,v and e keep their initial values
384	!-- (ug(nzt+1), vg(nzt+1), 0), at the bottom boundary the mirror
385	!-- boundary condition applies to u and v.
386	!-- The boundary condition for e is set further below ( (u/cm)*2 ).
387	u1d_p(nzb) = -u1d_p(nzb+1)
388	v1d_p(nzb) = -v1d_p(nzb+1)
389
390	!
391	!-- If necessary, apply the time filter
392	IF ( asselin_filter_factor /= 0.0 .AND. &
393	timestep_scheme(1:5) /= 'runge' ) THEN
394
395	u1d = u1d + asselin_filter_factor * ( u1d_p - 2.0 * u1d + u1d_m )
396	v1d = v1d + asselin_filter_factor * ( v1d_p - 2.0 * v1d + v1d_m )
397
398	IF ( .NOT. constant_diffusion ) THEN
399	e1d = e1d + asselin_filter_factor * &
400	( e1d_p - 2.0 * e1d + e1d_m )
401	ENDIF
402
403	ENDIF
404
405	!
406	!-- Swap the time levels in preparation for the next time step.
407	IF ( timestep_scheme(1:4) == 'leap' ) THEN
408	u1d_m = u1d
409	v1d_m = v1d
410	IF ( .NOT. constant_diffusion ) THEN
411	e1d_m = e1d
412	kh1d_m = kh1d ! The old diffusion quantities are required for
413	km1d_m = km1d ! explicit diffusion in the leap-frog scheme.
414	l1d_m = l1d
415	IF ( prandtl_layer ) THEN
416	usws1d_m = usws1d
417	vsws1d_m = vsws1d
418	ENDIF
419	ENDIF
420	ENDIF
421	u1d = u1d_p
422	v1d = v1d_p
423	IF ( .NOT. constant_diffusion ) THEN
424	e1d = e1d_p
425	ENDIF
426
427	!
428	!-- Compute diffusion quantities
429	IF ( .NOT. constant_diffusion ) THEN
430
431	!
432	!-- First compute the vertical fluxes in the Prandtl-layer
433	IF ( prandtl_layer ) THEN
434	!
435	!-- Compute theta* using Rif numbers of the previous time step
436	IF ( rif1d(1) >= 0.0 ) THEN
437	!
438	!-- Stable stratification
439	ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) / &
440	( LOG( zu(nzb+1) / z01d ) + 5.0 * rif1d(nzb+1) * &
441	( zu(nzb+1) - z01d ) / zu(nzb+1) &
442	)
443	ELSE
444	!
445	!-- Unstable stratification
446	a = SQRT( 1.0 - 16.0 * rif1d(nzb+1) )
447	b = SQRT( 1.0 - 16.0 * rif1d(nzb+1) / zu(nzb+1) * z01d )
448	!
449	!-- In the borderline case the formula for stable stratification
450	!-- must be applied, because otherwise a zero division would
451	!-- occur in the argument of the logarithm.
452	IF ( a == 0.0 .OR. b == 0.0 ) THEN
453	ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) / &
454	( LOG( zu(nzb+1) / z01d ) + 5.0 * rif1d(nzb+1) * &
455	( zu(nzb+1) - z01d ) / zu(nzb+1) &
456	)
457	ELSE
458	ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) / &
459	LOG( (a-1.0) / (a+1.0) * (b+1.0) / (b-1.0) )
460	ENDIF
461	ENDIF
462
463	ENDIF ! prandtl_layer
464
465	!
466	!-- Compute the Richardson-flux numbers,
467	!-- first at the top of the Prandtl-layer using u* of the previous
468	!-- time step (+1E-30, if u* = 0), then in the remaining area. There
469	!-- the rif-numbers of the previous time step are used.
470
471	IF ( prandtl_layer ) THEN
472	IF ( .NOT. humidity ) THEN
473	pt_0 = pt_init(nzb+1)
474	flux = ts1d
475	ELSE
476	pt_0 = pt_init(nzb+1) * ( 1.0 + 0.61 * q_init(nzb+1) )
477	flux = ts1d + 0.61 * pt_init(k) * qs1d
478	ENDIF
479	rif1d(nzb+1) = zu(nzb+1) * kappa * g * flux / &
480	( pt_0 * ( us1d**2 + 1E-30 ) )
481	ENDIF
482
483	DO k = nzb_diff, nzt
484	IF ( .NOT. humidity ) THEN
485	pt_0 = pt_init(k)
486	flux = ( pt_init(k+1) - pt_init(k-1) ) * dd2zu(k)
487	ELSE
488	pt_0 = pt_init(k) * ( 1.0 + 0.61 * q_init(k) )
489	flux = ( ( pt_init(k+1) - pt_init(k-1) ) &
490	+ 0.61 * pt_init(k) * ( q_init(k+1) - q_init(k-1) )&
491	) * dd2zu(k)
492	ENDIF
493	IF ( rif1d(k) >= 0.0 ) THEN
494	rif1d(k) = g / pt_0 * flux / &
495	( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2 &
496	+ ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2 &
497	+ 1E-30 &
498	)
499	ELSE
500	rif1d(k) = g / pt_0 * flux / &
501	( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2 &
502	+ ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2 &
503	+ 1E-30 &
504	) * ( 1.0 - 16.0 * rif1d(k) )**0.25
505	ENDIF
506	ENDDO
507	!
508	!-- Richardson-numbers must remain restricted to a realistic value
509	!-- range. It is exceeded excessively for very small velocities
510	!-- (u,v --> 0).
511	WHERE ( rif1d < rif_min ) rif1d = rif_min
512	WHERE ( rif1d > rif_max ) rif1d = rif_max
513
514	!
515	!-- Compute u* from the absolute velocity value
516	IF ( prandtl_layer ) THEN
517	uv_total = SQRT( u1d(nzb+1)2 + v1d(nzb+1)2 )
518
519	IF ( rif1d(nzb+1) >= 0.0 ) THEN
520	!
521	!-- Stable stratification
522	us1d = kappa * uv_total / ( &
523	LOG( zu(nzb+1) / z01d ) + 5.0 * rif1d(nzb+1) * &
524	( zu(nzb+1) - z01d ) / zu(nzb+1) &
525	)
526	ELSE
527	!
528	!-- Unstable stratification
529	a = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif1d(nzb+1) ) )
530	b = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif1d(nzb+1) / zu(nzb+1) &
531	* z01d ) )
532	!
533	!-- In the borderline case the formula for stable stratification
534	!-- must be applied, because otherwise a zero division would
535	!-- occur in the argument of the logarithm.
536	IF ( a == 1.0 .OR. b == 1.0 ) THEN
537	us1d = kappa * uv_total / ( &
538	LOG( zu(nzb+1) / z01d ) + &
539	5.0 * rif1d(nzb+1) * ( zu(nzb+1) - z01d ) / &
540	zu(nzb+1) )
541	ELSE
542	us1d = kappa * uv_total / ( &
543	LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) +&
544	2.0 * ( ATAN( b ) - ATAN( a ) ) &
545	)
546	ENDIF
547	ENDIF
548
549	!
550	!-- Compute the momentum fluxes for the diffusion terms
551	usws1d = - u1d(nzb+1) / uv_total * us1d**2
552	vsws1d = - v1d(nzb+1) / uv_total * us1d**2
553
554	!
555	!-- Boundary condition for the turbulent kinetic energy at the top
556	!-- of the Prandtl-layer. c_m = 0.4 according to Detering.
557	!-- Additional Neumann condition de/dz = 0 at nzb is set to ensure
558	!-- compatibility with the 3D model.
559	IF ( ibc_e_b == 2 ) THEN
560	e1d(nzb+1) = ( us1d / 0.1 )**2
561	! e1d(nzb+1) = ( us1d / 0.4 )**2 !not used so far, see also
562	!prandtl_fluxes
563	ENDIF
564	e1d(nzb) = e1d(nzb+1)
565
566	IF ( humidity .OR. passive_scalar ) THEN
567	!
568	!-- Compute q*
569	IF ( rif1d(1) >= 0.0 ) THEN
570	!
571	!-- Stable stratification
572	qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) / &
573	( LOG( zu(nzb+1) / z01d ) + 5.0 * rif1d(nzb+1) * &
574	( zu(nzb+1) - z01d ) / zu(nzb+1) &
575	)
576	ELSE
577	!
578	!-- Unstable stratification
579	a = SQRT( 1.0 - 16.0 * rif1d(nzb+1) )
580	b = SQRT( 1.0 - 16.0 * rif1d(nzb+1) / zu(nzb+1) * z01d )
581	!
582	!-- In the borderline case the formula for stable stratification
583	!-- must be applied, because otherwise a zero division would
584	!-- occur in the argument of the logarithm.
585	IF ( a == 1.0 .OR. b == 1.0 ) THEN
586	qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) / &
587	( LOG( zu(nzb+1) / z01d ) + 5.0 * rif1d(nzb+1) * &
588	( zu(nzb+1) - z01d ) / zu(nzb+1) &
589	)
590	ELSE
591	qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) / &
592	LOG( (a-1.0) / (a+1.0) * (b+1.0) / (b-1.0) )
593	ENDIF
594	ENDIF
595	ELSE
596	qs1d = 0.0
597	ENDIF
598
599	ENDIF ! prandtl_layer
600
601	!
602	!-- Compute the diabatic mixing length
603	IF ( mixing_length_1d == 'blackadar' ) THEN
604	DO k = nzb+1, nzt
605	IF ( rif1d(k) >= 0.0 ) THEN
606	l1d(k) = l_black(k) / ( 1.0 + 5.0 * rif1d(k) )
607	ELSE
608	l1d(k) = l_black(k) * ( 1.0 - 16.0 * rif1d(k) )**0.25
609	ENDIF
610	l1d(k) = l_black(k)
611	ENDDO
612
613	ELSEIF ( mixing_length_1d == 'as_in_3d_model' ) THEN
614	DO k = nzb+1, nzt
615	dpt_dz = ( pt_init(k+1) - pt_init(k-1) ) * dd2zu(k)
616	IF ( dpt_dz > 0.0 ) THEN
617	l_stable = 0.76 * SQRT( e1d(k) ) / &
618	SQRT( g / pt_init(k) * dpt_dz ) + 1E-5
619	ELSE
620	l_stable = l_grid(k)
621	ENDIF
622	l1d(k) = MIN( l_grid(k), l_stable )
623	ENDDO
624	ENDIF
625
626	!
627	!-- Adjust mixing length to the prandtl mixing length
628	IF ( adjust_mixing_length .AND. prandtl_layer ) THEN
629	k = nzb+1
630	IF ( rif1d(k) >= 0.0 ) THEN
631	l1d(k) = MIN( l1d(k), kappa * zu(k) / ( 1.0 + 5.0 * &
632	rif1d(k) ) )
633	ELSE
634	l1d(k) = MIN( l1d(k), kappa * zu(k) * &
635	SQRT( SQRT( 1.0 - 16.0 * rif1d(k) ) ) )
636	ENDIF
637	ENDIF
638
639	!
640	!-- Compute the diffusion coefficients for momentum via the
641	!-- corresponding Prandtl-layer relationship and according to
642	!-- Prandtl-Kolmogorov, respectively. The unstable stratification is
643	!-- computed via the adiabatic mixing length, for the unstability has
644	!-- already been taken account of via the TKE (cf. also Diss.).
645	IF ( prandtl_layer ) THEN
646	IF ( rif1d(nzb+1) >= 0.0 ) THEN
647	km1d(nzb+1) = us1d * kappa * zu(nzb+1) / &
648	( 1.0 + 5.0 * rif1d(nzb+1) )
649	ELSE
650	km1d(nzb+1) = us1d * kappa * zu(nzb+1) * &
651	( 1.0 - 16.0 * rif1d(nzb+1) )**0.25
652	ENDIF
653	ENDIF
654	DO k = nzb_diff, nzt
655	! km1d(k) = 0.4 * SQRT( e1d(k) ) !changed: adjustment to 3D-model
656	km1d(k) = 0.1 * SQRT( e1d(k) )
657	IF ( rif1d(k) >= 0.0 ) THEN
658	km1d(k) = km1d(k) * l1d(k)
659	ELSE
660	km1d(k) = km1d(k) * l_black(k)
661	ENDIF
662	ENDDO
663
664	!
665	!-- Add damping layer
666	DO k = damp_level_ind_1d+1, nzt+1
667	km1d(k) = 1.1 * km1d(k-1)
668	km1d(k) = MIN( km1d(k), 10.0 )
669	ENDDO
670
671	!
672	!-- Compute the diffusion coefficient for heat via the relationship
673	!-- kh = phim / phih * km
674	DO k = nzb+1, nzt
675	IF ( rif1d(k) >= 0.0 ) THEN
676	kh1d(k) = km1d(k)
677	ELSE
678	kh1d(k) = km1d(k) * ( 1.0 - 16.0 * rif1d(k) )**0.25
679	ENDIF
680	ENDDO
681
682	ENDIF ! .NOT. constant_diffusion
683
684	!
685	!-- The Runge-Kutta scheme needs the recent diffusion quantities
686	IF ( timestep_scheme(1:5) == 'runge' ) THEN
687	u1d_m = u1d
688	v1d_m = v1d
689	IF ( .NOT. constant_diffusion ) THEN
690	e1d_m = e1d
691	kh1d_m = kh1d
692	km1d_m = km1d
693	l1d_m = l1d
694	IF ( prandtl_layer ) THEN
695	usws1d_m = usws1d
696	vsws1d_m = vsws1d
697	ENDIF
698	ENDIF
699	ENDIF
700
701
702	ENDDO ! intermediate step loop
703
704	!
705	!-- Increment simulated time and output times
706	current_timestep_number_1d = current_timestep_number_1d + 1
707	simulated_time_1d = simulated_time_1d + dt_1d
708	simulated_time_chr = time_to_string( simulated_time_1d )
709	time_pr_1d = time_pr_1d + dt_1d
710	time_run_control_1d = time_run_control_1d + dt_1d
711
712	!
713	!-- Determine and print out quantities for run control
714	IF ( time_run_control_1d >= dt_run_control_1d ) THEN
715	CALL run_control_1d
716	time_run_control_1d = time_run_control_1d - dt_run_control_1d
717	ENDIF
718
719	!
720	!-- Profile output on file
721	IF ( time_pr_1d >= dt_pr_1d ) THEN
722	CALL print_1d_model
723	time_pr_1d = time_pr_1d - dt_pr_1d
724	ENDIF
725
726	!
727	!-- Determine size of next time step
728	CALL timestep_1d
729
730	ENDDO ! time loop
731
732
733	END SUBROUTINE time_integration_1d
734
735
736	SUBROUTINE run_control_1d
737
738	!------------------------------------------------------------------------------!
739	! Description:
740	! ------------
741	! Compute and print out quantities for run control of the 1D model.
742	!------------------------------------------------------------------------------!
743
744	USE constants
745	USE indices
746	USE model_1d
747	USE pegrid
748	USE control_parameters
749
750	IMPLICIT NONE
751
752	INTEGER :: k
753	REAL :: alpha, energy, umax, uv_total, vmax
754
755	!
756	!-- Output
757	IF ( myid == 0 ) THEN
758	!
759	!-- If necessary, write header
760	IF ( .NOT. run_control_header_1d ) THEN
761	CALL check_open( 15 )
762	WRITE ( 15, 100 )
763	run_control_header_1d = .TRUE.
764	ENDIF
765
766	!
767	!-- Compute control quantities
768	!-- grid level nzp is excluded due to mirror boundary condition
769	umax = 0.0; vmax = 0.0; energy = 0.0
770	DO k = nzb+1, nzt+1
771	umax = MAX( ABS( umax ), ABS( u1d(k) ) )
772	vmax = MAX( ABS( vmax ), ABS( v1d(k) ) )
773	energy = energy + 0.5 * ( u1d(k)2 + v1d(k)2 )
774	ENDDO
775	energy = energy / REAL( nzt - nzb + 1 )
776
777	uv_total = SQRT( u1d(nzb+1)2 + v1d(nzb+1)2 )
778	IF ( ABS( v1d(nzb+1) ) .LT. 1.0E-5 ) THEN
779	alpha = ACOS( SIGN( 1.0 , u1d(nzb+1) ) )
780	ELSE
781	alpha = ACOS( u1d(nzb+1) / uv_total )
782	IF ( v1d(nzb+1) <= 0.0 ) alpha = 2.0 * pi - alpha
783	ENDIF
784	alpha = alpha / ( 2.0 * pi ) * 360.0
785
786	WRITE ( 15, 101 ) current_timestep_number_1d, simulated_time_chr, &
787	dt_1d, umax, vmax, us1d, alpha, energy
788	!
789	!-- Write buffer contents to disc immediately
790	CALL local_flush( 15 )
791
792	ENDIF
793
794	!
795	!-- formats
796	100 FORMAT (///'1D-Zeitschrittkontrollausgaben:'/ &
797	&'------------------------------'// &
798	&'ITER. HH:MM:SS DT UMAX VMAX U* ALPHA ENERG.'/ &
799	&'-------------------------------------------------------------')
800	101 FORMAT (I5,2X,A9,1X,F6.2,2X,F6.2,1X,F6.2,2X,F5.3,2X,F5.1,2X,F7.2)
801
802
803	END SUBROUTINE run_control_1d
804
805
806
807	SUBROUTINE timestep_1d
808
809	!------------------------------------------------------------------------------!
810	! Description:
811	! ------------
812	! Compute the time step w.r.t. the diffusion criterion
813	!------------------------------------------------------------------------------!
814
815	USE arrays_3d
816	USE indices
817	USE model_1d
818	USE pegrid
819	USE control_parameters
820
821	IMPLICIT NONE
822
823	INTEGER :: k
824	REAL :: div, dt_diff, fac, percent_change, value
825
826
827	!
828	!-- Compute the currently feasible time step according to the diffusion
829	!-- criterion. At nzb+1 the half grid length is used.
830	IF ( timestep_scheme(1:4) == 'leap' ) THEN
831	fac = 0.25
832	ELSE
833	fac = 0.35
834	ENDIF
835	dt_diff = dt_max_1d
836	DO k = nzb+2, nzt
837	value = fac * dzu(k) * dzu(k) / ( km1d(k) + 1E-20 )
838	dt_diff = MIN( value, dt_diff )
839	ENDDO
840	value = fac * zu(nzb+1) * zu(nzb+1) / ( km1d(nzb+1) + 1E-20 )
841	dt_1d = MIN( value, dt_diff )
842
843	!
844	!-- Set flag when the time step becomes too small
845	IF ( dt_1d < ( 0.00001 * dt_max_1d ) ) THEN
846	stop_dt_1d = .TRUE.
847
848	WRITE( message_string, * ) 'timestep has exceeded the lower limit &', &
849	'dt_1d = ',dt_1d,' s simulation stopped!'
850	CALL message( 'timestep_1d', 'PA0192', 1, 2, 0, 6, 0 )
851
852	ENDIF
853
854	IF ( timestep_scheme(1:4) == 'leap' ) THEN
855
856	!
857	!-- The current time step will only be changed if the new time step exceeds
858	!-- its previous value by 5 % or falls 2 % below. After a time step
859	!-- reduction at least 30 iterations must be done with this value before a
860	!-- new reduction will be allowed again.
861	!-- The control parameters for application of Euler- or leap-frog schemes are
862	!-- set accordingly.
863	percent_change = dt_1d / old_dt_1d - 1.0
864	IF ( percent_change > 0.05 .OR. percent_change < -0.02 ) THEN
865
866	!
867	!-- Each time step increase is by at most 2 %
868	IF ( percent_change > 0.0 .AND. simulated_time_1d /= 0.0 ) THEN
869	dt_1d = 1.02 * old_dt_1d
870	ENDIF
871
872	!
873	!-- A more or less simple new time step value is obtained taking only the
874	!-- first two significant digits
875	div = 1000.0
876	DO WHILE ( dt_1d < div )
877	div = div / 10.0
878	ENDDO
879	dt_1d = NINT( dt_1d * 100.0 / div ) * div / 100.0
880
881	!
882	!-- Now the time step can be changed.
883	IF ( percent_change < 0.0 ) THEN
884	!
885	!-- Time step reduction
886	old_dt_1d = dt_1d
887	last_dt_change_1d = current_timestep_number_1d
888	ELSE
889	!
890	!-- Time step will only be increased if at least 30 iterations have
891	!-- been done since the previous time step change, and of course at
892	!-- simulation start, respectively.
893	IF ( current_timestep_number_1d >= last_dt_change_1d + 30 .OR. &
894	simulated_time_1d == 0.0 ) THEN
895	old_dt_1d = dt_1d
896	last_dt_change_1d = current_timestep_number_1d
897	ELSE
898	dt_1d = old_dt_1d
899	ENDIF
900	ENDIF
901	ELSE
902	!
903	!-- No time step change since the difference is too small
904	dt_1d = old_dt_1d
905	ENDIF
906
907	ELSE ! Runge-Kutta
908
909	!-- A more or less simple new time step value is obtained taking only the
910	!-- first two significant digits
911	div = 1000.0
912	DO WHILE ( dt_1d < div )
913	div = div / 10.0
914	ENDDO
915	dt_1d = NINT( dt_1d * 100.0 / div ) * div / 100.0
916
917	old_dt_1d = dt_1d
918	last_dt_change_1d = current_timestep_number_1d
919
920	ENDIF
921
922	END SUBROUTINE timestep_1d
923
924
925
926	SUBROUTINE print_1d_model
927
928	!------------------------------------------------------------------------------!
929	! Description:
930	! ------------
931	! List output of profiles from the 1D-model
932	!------------------------------------------------------------------------------!
933
934	USE arrays_3d
935	USE indices
936	USE model_1d
937	USE pegrid
938	USE control_parameters
939
940	IMPLICIT NONE
941
942
943	INTEGER :: k
944
945
946	IF ( myid == 0 ) THEN
947	!
948	!-- Open list output file for profiles from the 1D-model
949	CALL check_open( 17 )
950
951	!
952	!-- Write Header
953	WRITE ( 17, 100 ) TRIM( run_description_header ), &
954	TRIM( simulated_time_chr )
955	WRITE ( 17, 101 )
956
957	!
958	!-- Write the values
959	WRITE ( 17, 102 )
960	WRITE ( 17, 101 )
961	DO k = nzt+1, nzb, -1
962	WRITE ( 17, 103) k, zu(k), u1d(k), v1d(k), pt_init(k), e1d(k), &
963	rif1d(k), km1d(k), kh1d(k), l1d(k), zu(k), k
964	ENDDO
965	WRITE ( 17, 101 )
966	WRITE ( 17, 102 )
967	WRITE ( 17, 101 )
968
969	!
970	!-- Write buffer contents to disc immediately
971	CALL local_flush( 17 )
972
973	ENDIF
974
975	!
976	!-- Formats
977	100 FORMAT (//1X,A/1X,10('-')/' 1d-model profiles'/ &
978	' Time: ',A)
979	101 FORMAT (1X,79('-'))
980	102 FORMAT (' k zu u v pt e rif Km Kh ', &
981	'l zu k')
982	103 FORMAT (1X,I4,1X,F7.1,1X,F6.2,1X,F6.2,1X,F6.2,1X,F6.2,1X,F5.2,1X,F5.2, &
983	1X,F5.2,1X,F6.2,1X,F7.1,2X,I4)
984
985
986	END SUBROUTINE print_1d_model

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

Context Navigation

source: palm/trunk/SOURCE/init_1d_model.f90 @ 419

Download in other formats: