Home

Context Navigation

init_1d_model.f90 @ 516

Last change on this file since 516 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

Rev	Line
[1]	1	SUBROUTINE init_1d_model
	2
	3	!------------------------------------------------------------------------------!
[254]	4	! Current revisions:
[1]	5	! -----------------
[392]	6	!
[1]	7	!
	8	! Former revisions:
	9	! -----------------
[3]	10	! $Id: init_1d_model.f90 392 2009-09-24 10:39:14Z raasch $
[77]	11	!
[392]	12	! 254 2009-03-05 15:33:42Z heinze
	13	! Output of messages replaced by message handling routine.
	14	!
[198]	15	! 184 2008-08-04 15:53:39Z letzel
	16	! provisional solution for run_control_1d output: add 'CALL check_open( 15 )'
	17	!
[139]	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	!
[83]	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	!
[77]	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	!
[3]	30	! RCS Log replace by Id keyword, revision history cleaned up
	31	!
[1]	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
[135]	88	lambda = 2.7E-4 * SQRT( ug(nzt+1)2 + vg(nzt+1)2 ) / &
	89	ABS( f ) + 1E-10
[1]	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
[75]	147	IF ( humidity .OR. passive_scalar ) qs1d = 0.0
[1]	148
	149	!
[46]	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	!
[1]	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) )
[75]	238	IF ( .NOT. humidity ) THEN
[1]	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) )
[75]	278	IF ( .NOT. humidity ) THEN
[1]	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
[75]	472	IF ( .NOT. humidity ) THEN
[1]	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
[75]	484	IF ( .NOT. humidity ) THEN
[1]	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
[75]	566	IF ( humidity .OR. passive_scalar ) THEN
[1]	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
[184]	761	CALL check_open( 15 )
[1]	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
[82]	790	CALL local_flush( 15 )
[1]	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.
[254]	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
[1]	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
[82]	971	CALL local_flush( 17 )
[1]	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 @ 516

Download in other formats: