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

Last change on this file since 392 was 392, checked in by raasch, 12 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 raasch $
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
796100 FORMAT (///'1D-Zeitschrittkontrollausgaben:'/ &
797              &'------------------------------'// &
798           &'ITER.  HH:MM:SS    DT      UMAX   VMAX    U*   ALPHA   ENERG.'/ &
799           &'-------------------------------------------------------------')
800101 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
977100 FORMAT (//1X,A/1X,10('-')/' 1d-model profiles'/ &
978            ' Time: ',A)
979101 FORMAT (1X,79('-'))
980102 FORMAT ('   k     zu      u      v     pt      e    rif    Km    Kh     ', &
981            'l      zu      k')
982103 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.