source: palm/trunk/SOURCE/init_3d_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: 55.0 KB
Line 
1#if defined( __ibmy_special )
2@PROCESS NOOPTimize
3#endif
4 SUBROUTINE init_3d_model
5
6!------------------------------------------------------------------------------!
7! Current revisions:
8! -----------------
9!
10!
11! Former revisions:
12! -----------------
13! $Id: init_3d_model.f90 392 2009-09-24 10:39:14Z raasch $
14!
15! 388 2009-09-23 09:40:33Z raasch
16! Initialization of prho added.
17! bugfix: correction of initial volume flow for non-flat topography
18! bugfix: zero initialization of arrays within buildings for 'cyclic_fill'
19! bugfix: avoid that ngp_2dh_s_inner becomes zero
20! initializing_actions='read_data_for_recycling' renamed to 'cyclic_fill', now
21! independent of turbulent_inflow
22! Output of messages replaced by message handling routine.
23! Set the starting level and the vertical smoothing factor used for
24! the external pressure gradient
25! +conserve_volume_flow_mode: 'default', 'initial_profiles', 'inflow_profile'
26! and 'bulk_velocity'
27! If the inversion height calculated by the prerun is zero,
28! inflow_damping_height must be explicitly specified.
29!
30! 181 2008-07-30 07:07:47Z raasch
31! bugfix: zero assignments to tendency arrays in case of restarts,
32! further extensions and modifications in the initialisation of the plant
33! canopy model,
34! allocation of hom_sum moved to parin, initialization of spectrum_x|y directly
35! after allocating theses arrays,
36! read data for recycling added as new initialization option,
37! dummy allocation for diss
38!
39! 138 2007-11-28 10:03:58Z letzel
40! New counter ngp_2dh_s_inner.
41! Allow new case bc_uv_t = 'dirichlet_0' for channel flow.
42! Corrected calculation of initial volume flow for 'set_1d-model_profiles' and
43! 'set_constant_profiles' in case of buildings in the reference cross-sections.
44!
45! 108 2007-08-24 15:10:38Z letzel
46! Flux initialization in case of coupled runs, +momentum fluxes at top boundary,
47! +arrays for phase speed c_u, c_v, c_w, indices for u|v|w_m_l|r changed
48! +qswst_remote in case of atmosphere model with humidity coupled to ocean
49! Rayleigh damping for ocean, optionally calculate km and kh from initial
50! TKE e_init
51!
52! 97 2007-06-21 08:23:15Z raasch
53! Initialization of salinity, call of init_ocean
54!
55! 87 2007-05-22 15:46:47Z raasch
56! var_hom and var_sum renamed pr_palm
57!
58! 75 2007-03-22 09:54:05Z raasch
59! Arrays for radiation boundary conditions are allocated (u_m_l, u_m_r, etc.),
60! bugfix for cases with the outflow damping layer extending over more than one
61! subdomain, moisture renamed humidity,
62! new initializing action "by_user" calls user_init_3d_model,
63! precipitation_amount/rate, ts_value are allocated, +module netcdf_control,
64! initial velocities at nzb+1 are regarded for volume
65! flow control in case they have been set zero before (to avoid small timesteps)
66! -uvmean_outflow, uxrp, vynp eliminated
67!
68! 19 2007-02-23 04:53:48Z raasch
69! +handling of top fluxes
70!
71! RCS Log replace by Id keyword, revision history cleaned up
72!
73! Revision 1.49  2006/08/22 15:59:07  raasch
74! No optimization of this file on the ibmy (Yonsei Univ.)
75!
76! Revision 1.1  1998/03/09 16:22:22  raasch
77! Initial revision
78!
79!
80! Description:
81! ------------
82! Allocation of arrays and initialization of the 3D model via
83! a) pre-run the 1D model
84! or
85! b) pre-set constant linear profiles
86! or
87! c) read values of a previous run
88!------------------------------------------------------------------------------!
89
90    USE arrays_3d
91    USE averaging
92    USE cloud_parameters
93    USE constants
94    USE control_parameters
95    USE cpulog
96    USE indices
97    USE interfaces
98    USE model_1d
99    USE netcdf_control
100    USE particle_attributes
101    USE pegrid
102    USE profil_parameter
103    USE random_function_mod
104    USE statistics
105
106    IMPLICIT NONE
107
108    INTEGER ::  i, ind_array(1), j, k, sr
109
110    INTEGER, DIMENSION(:), ALLOCATABLE ::  ngp_2dh_l, ngp_3d_inner_l
111
112    INTEGER, DIMENSION(:,:), ALLOCATABLE ::  ngp_2dh_outer_l,  &
113         ngp_2dh_s_inner_l
114
115    REAL ::  a, b
116
117    REAL, DIMENSION(1:2) ::  volume_flow_area_l, volume_flow_initial_l
118
119
120!
121!-- Allocate arrays
122    ALLOCATE( ngp_2dh(0:statistic_regions), ngp_2dh_l(0:statistic_regions), &
123              ngp_3d(0:statistic_regions),                                  &
124              ngp_3d_inner(0:statistic_regions),                            &
125              ngp_3d_inner_l(0:statistic_regions),                          &
126              sums_divnew_l(0:statistic_regions),                           &
127              sums_divold_l(0:statistic_regions) )
128    ALLOCATE( dp_smooth_factor(nzb:nzt), rdf(nzb+1:nzt) )
129    ALLOCATE( ngp_2dh_outer(nzb:nzt+1,0:statistic_regions),                 &
130              ngp_2dh_outer_l(nzb:nzt+1,0:statistic_regions),               &
131              ngp_2dh_s_inner(nzb:nzt+1,0:statistic_regions),               &
132              ngp_2dh_s_inner_l(nzb:nzt+1,0:statistic_regions),             &
133              rmask(nys-1:nyn+1,nxl-1:nxr+1,0:statistic_regions),           &
134              sums(nzb:nzt+1,pr_palm+max_pr_user),                          &
135              sums_l(nzb:nzt+1,pr_palm+max_pr_user,0:threads_per_task-1),   &
136              sums_l_l(nzb:nzt+1,0:statistic_regions,0:threads_per_task-1), &
137              sums_up_fraction_l(10,3,0:statistic_regions),                 &
138              sums_wsts_bc_l(nzb:nzt+1,0:statistic_regions),                &
139              ts_value(var_ts,0:statistic_regions) )
140    ALLOCATE( km_damp_x(nxl-1:nxr+1), km_damp_y(nys-1:nyn+1) )
141
142    ALLOCATE( rif_1(nys-1:nyn+1,nxl-1:nxr+1), shf_1(nys-1:nyn+1,nxl-1:nxr+1), &
143              ts(nys-1:nyn+1,nxl-1:nxr+1), tswst_1(nys-1:nyn+1,nxl-1:nxr+1),  &
144              us(nys-1:nyn+1,nxl-1:nxr+1), usws_1(nys-1:nyn+1,nxl-1:nxr+1),   &
145              uswst_1(nys-1:nyn+1,nxl-1:nxr+1),                               &
146              vsws_1(nys-1:nyn+1,nxl-1:nxr+1),                                &
147              vswst_1(nys-1:nyn+1,nxl-1:nxr+1), z0(nys-1:nyn+1,nxl-1:nxr+1) )
148
149    IF ( timestep_scheme(1:5) /= 'runge' )  THEN
150!
151!--    Leapfrog scheme needs two timelevels of diffusion quantities
152       ALLOCATE( rif_2(nys-1:nyn+1,nxl-1:nxr+1),   &
153                 shf_2(nys-1:nyn+1,nxl-1:nxr+1),   &
154                 tswst_2(nys-1:nyn+1,nxl-1:nxr+1), &
155                 usws_2(nys-1:nyn+1,nxl-1:nxr+1),  &
156                 uswst_2(nys-1:nyn+1,nxl-1:nxr+1), &
157                 vswst_2(nys-1:nyn+1,nxl-1:nxr+1), &
158                 vsws_2(nys-1:nyn+1,nxl-1:nxr+1) )
159    ENDIF
160
161    ALLOCATE( d(nzb+1:nzta,nys:nyna,nxl:nxra),         &
162              e_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
163              e_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
164              e_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
165              kh_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
166              km_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
167              p(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),    &
168              pt_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
169              pt_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
170              pt_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
171              tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
172              u_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
173              u_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
174              u_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
175              v_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
176              v_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
177              v_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
178              w_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
179              w_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
180              w_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
181
182    IF ( timestep_scheme(1:5) /= 'runge' )  THEN
183       ALLOCATE( kh_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
184                 km_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
185    ENDIF
186
187    IF ( humidity  .OR.  passive_scalar ) THEN
188!
189!--    2D-humidity/scalar arrays
190       ALLOCATE ( qs(nys-1:nyn+1,nxl-1:nxr+1),     &
191                  qsws_1(nys-1:nyn+1,nxl-1:nxr+1), &
192                  qswst_1(nys-1:nyn+1,nxl-1:nxr+1) )
193
194       IF ( timestep_scheme(1:5) /= 'runge' )  THEN
195          ALLOCATE( qsws_2(nys-1:nyn+1,nxl-1:nxr+1), &
196                    qswst_2(nys-1:nyn+1,nxl-1:nxr+1) )
197       ENDIF
198!
199!--    3D-humidity/scalar arrays
200       ALLOCATE( q_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
201                 q_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
202                 q_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
203
204!
205!--    3D-arrays needed for humidity only
206       IF ( humidity )  THEN
207          ALLOCATE( vpt_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
208
209          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
210             ALLOCATE( vpt_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
211          ENDIF
212
213          IF ( cloud_physics ) THEN
214!
215!--          Liquid water content
216             ALLOCATE ( ql_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
217!
218!--          Precipitation amount and rate (only needed if output is switched)
219             ALLOCATE( precipitation_amount(nys-1:nyn+1,nxl-1:nxr+1), &
220                       precipitation_rate(nys-1:nyn+1,nxl-1:nxr+1) )
221          ENDIF
222
223          IF ( cloud_droplets )  THEN
224!
225!--          Liquid water content, change in liquid water content,
226!--          real volume of particles (with weighting), volume of particles
227             ALLOCATE ( ql_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
228                        ql_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
229                        ql_v(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
230                        ql_vp(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
231          ENDIF
232
233       ENDIF
234
235    ENDIF
236
237    IF ( ocean )  THEN
238       ALLOCATE( saswsb_1(nys-1:nyn+1,nxl-1:nxr+1), &
239                 saswst_1(nys-1:nyn+1,nxl-1:nxr+1) )
240       ALLOCATE( prho_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
241                 rho_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
242                 sa_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),   &
243                 sa_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),   &
244                 sa_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
245       prho => prho_1
246       rho  => rho_1  ! routines calc_mean_profile and diffusion_e require
247                      ! density to be apointer
248       IF ( humidity_remote )  THEN
249          ALLOCATE( qswst_remote(nys-1:nyn+1,nxl-1:nxr+1) )
250          qswst_remote = 0.0
251       ENDIF
252    ENDIF
253
254!
255!-- 3D-array for storing the dissipation, needed for calculating the sgs
256!-- particle velocities
257    IF ( use_sgs_for_particles )  THEN
258       ALLOCATE ( diss(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
259    ELSE
260       ALLOCATE ( diss(2,2,2) )  ! required because diss is used as a
261                                 ! formal parameter
262    ENDIF
263
264    IF ( dt_dosp /= 9999999.9 )  THEN
265       ALLOCATE( spectrum_x( 1:nx/2, 1:10, 1:10 ), &
266                 spectrum_y( 1:ny/2, 1:10, 1:10 ) )
267       spectrum_x = 0.0
268       spectrum_y = 0.0
269    ENDIF
270
271!
272!-- 3D-arrays for the leaf area density and the canopy drag coefficient
273    IF ( plant_canopy ) THEN
274       ALLOCATE ( lad_s(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
275                  lad_u(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
276                  lad_v(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
277                  lad_w(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
278                  cdc(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
279
280       IF ( passive_scalar ) THEN
281          ALLOCATE ( sls(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),   &
282                     sec(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) ) 
283       ENDIF
284
285       IF ( cthf /= 0.0 ) THEN
286          ALLOCATE ( lai(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),   &
287                     canopy_heat_flux(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
288       ENDIF
289
290    ENDIF
291
292!
293!-- 4D-array for storing the Rif-values at vertical walls
294    IF ( topography /= 'flat' )  THEN
295       ALLOCATE( rif_wall(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1,1:4) )
296       rif_wall = 0.0
297    ENDIF
298
299!
300!-- Velocities at nzb+1 needed for volume flow control
301    IF ( conserve_volume_flow )  THEN
302       ALLOCATE( u_nzb_p1_for_vfc(nys:nyn), v_nzb_p1_for_vfc(nxl:nxr) )
303       u_nzb_p1_for_vfc = 0.0
304       v_nzb_p1_for_vfc = 0.0
305    ENDIF
306
307!
308!-- Arrays to store velocity data from t-dt and the phase speeds which
309!-- are needed for radiation boundary conditions
310    IF ( outflow_l )  THEN
311       ALLOCATE( u_m_l(nzb:nzt+1,nys-1:nyn+1,1:2), &
312                 v_m_l(nzb:nzt+1,nys-1:nyn+1,0:1), &
313                 w_m_l(nzb:nzt+1,nys-1:nyn+1,0:1) )
314    ENDIF
315    IF ( outflow_r )  THEN
316       ALLOCATE( u_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx), &
317                 v_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx), &
318                 w_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx) )
319    ENDIF
320    IF ( outflow_l  .OR.  outflow_r )  THEN
321       ALLOCATE( c_u(nzb:nzt+1,nys-1:nyn+1), c_v(nzb:nzt+1,nys-1:nyn+1), &
322                 c_w(nzb:nzt+1,nys-1:nyn+1) )
323    ENDIF
324    IF ( outflow_s )  THEN
325       ALLOCATE( u_m_s(nzb:nzt+1,0:1,nxl-1:nxr+1), &
326                 v_m_s(nzb:nzt+1,1:2,nxl-1:nxr+1), &
327                 w_m_s(nzb:nzt+1,0:1,nxl-1:nxr+1) )
328    ENDIF
329    IF ( outflow_n )  THEN
330       ALLOCATE( u_m_n(nzb:nzt+1,ny-1:ny,nxl-1:nxr+1), &
331                 v_m_n(nzb:nzt+1,ny-1:ny,nxl-1:nxr+1), &
332                 w_m_n(nzb:nzt+1,ny-1:ny,nxl-1:nxr+1) )
333    ENDIF
334    IF ( outflow_s  .OR.  outflow_n )  THEN
335       ALLOCATE( c_u(nzb:nzt+1,nxl-1:nxr+1), c_v(nzb:nzt+1,nxl-1:nxr+1), &
336                 c_w(nzb:nzt+1,nxl-1:nxr+1) )
337    ENDIF
338
339!
340!-- Initial assignment of the pointers
341    IF ( timestep_scheme(1:5) /= 'runge' )  THEN
342
343       rif_m   => rif_1;    rif   => rif_2
344       shf_m   => shf_1;    shf   => shf_2
345       tswst_m => tswst_1;  tswst => tswst_2
346       usws_m  => usws_1;   usws  => usws_2
347       uswst_m => uswst_1;  uswst => uswst_2
348       vsws_m  => vsws_1;   vsws  => vsws_2
349       vswst_m => vswst_1;  vswst => vswst_2
350       e_m  => e_1;   e  => e_2;   e_p  => e_3;   te_m  => e_3
351       kh_m => kh_1;  kh => kh_2
352       km_m => km_1;  km => km_2
353       pt_m => pt_1;  pt => pt_2;  pt_p => pt_3;  tpt_m => pt_3
354       u_m  => u_1;   u  => u_2;   u_p  => u_3;   tu_m  => u_3
355       v_m  => v_1;   v  => v_2;   v_p  => v_3;   tv_m  => v_3
356       w_m  => w_1;   w  => w_2;   w_p  => w_3;   tw_m  => w_3
357
358       IF ( humidity  .OR.  passive_scalar )  THEN
359          qsws_m  => qsws_1;   qsws  => qsws_2
360          qswst_m => qswst_1;  qswst => qswst_2
361          q_m    => q_1;     q    => q_2;     q_p => q_3;     tq_m => q_3
362          IF ( humidity )        vpt_m  => vpt_1;   vpt  => vpt_2
363          IF ( cloud_physics )   ql   => ql_1
364          IF ( cloud_droplets )  THEN
365             ql   => ql_1
366             ql_c => ql_2
367          ENDIF
368       ENDIF
369
370    ELSE
371
372       rif   => rif_1
373       shf   => shf_1
374       tswst => tswst_1
375       usws  => usws_1
376       uswst => uswst_1
377       vsws  => vsws_1
378       vswst => vswst_1
379       e     => e_1;   e_p  => e_2;   te_m  => e_3;   e_m  => e_3
380       kh    => kh_1
381       km    => km_1
382       pt    => pt_1;  pt_p => pt_2;  tpt_m => pt_3;  pt_m => pt_3
383       u     => u_1;   u_p  => u_2;   tu_m  => u_3;   u_m  => u_3
384       v     => v_1;   v_p  => v_2;   tv_m  => v_3;   v_m  => v_3
385       w     => w_1;   w_p  => w_2;   tw_m  => w_3;   w_m  => w_3
386
387       IF ( humidity  .OR.  passive_scalar )  THEN
388          qsws   => qsws_1
389          qswst  => qswst_1
390          q      => q_1;     q_p  => q_2;     tq_m  => q_3;    q_m => q_3
391          IF ( humidity )        vpt  => vpt_1
392          IF ( cloud_physics )   ql   => ql_1
393          IF ( cloud_droplets )  THEN
394             ql   => ql_1
395             ql_c => ql_2
396          ENDIF
397       ENDIF
398
399       IF ( ocean )  THEN
400          saswsb => saswsb_1
401          saswst => saswst_1
402          sa     => sa_1;    sa_p => sa_2;    tsa_m => sa_3
403       ENDIF
404
405    ENDIF
406
407!
408!-- Initialize model variables
409    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.  &
410         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
411!
412!--    First model run of a possible job queue.
413!--    Initial profiles of the variables must be computes.
414       IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
415!
416!--       Use solutions of the 1D model as initial profiles,
417!--       start 1D model
418          CALL init_1d_model
419!
420!--       Transfer initial profiles to the arrays of the 3D model
421          DO  i = nxl-1, nxr+1
422             DO  j = nys-1, nyn+1
423                e(:,j,i)  = e1d
424                kh(:,j,i) = kh1d
425                km(:,j,i) = km1d
426                pt(:,j,i) = pt_init
427                u(:,j,i)  = u1d
428                v(:,j,i)  = v1d
429             ENDDO
430          ENDDO
431
432          IF ( humidity  .OR.  passive_scalar )  THEN
433             DO  i = nxl-1, nxr+1
434                DO  j = nys-1, nyn+1
435                   q(:,j,i) = q_init
436                ENDDO
437             ENDDO
438          ENDIF
439
440          IF ( .NOT. constant_diffusion )  THEN
441             DO  i = nxl-1, nxr+1
442                DO  j = nys-1, nyn+1
443                   e(:,j,i)  = e1d
444                ENDDO
445             ENDDO
446!
447!--          Store initial profiles for output purposes etc.
448             hom(:,1,25,:) = SPREAD( l1d, 2, statistic_regions+1 )
449
450             IF ( prandtl_layer )  THEN
451                rif  = rif1d(nzb+1)
452                ts   = 0.0  ! could actually be computed more accurately in the
453                            ! 1D model. Update when opportunity arises.
454                us   = us1d
455                usws = usws1d
456                vsws = vsws1d
457             ELSE
458                ts   = 0.0  ! must be set, because used in
459                rif  = 0.0  ! flowste
460                us   = 0.0
461                usws = 0.0
462                vsws = 0.0
463             ENDIF
464
465          ELSE
466             e    = 0.0  ! must be set, because used in
467             rif  = 0.0  ! flowste
468             ts   = 0.0
469             us   = 0.0
470             usws = 0.0
471             vsws = 0.0
472          ENDIF
473          uswst = top_momentumflux_u
474          vswst = top_momentumflux_v
475
476!
477!--       In every case qs = 0.0 (see also pt)
478!--       This could actually be computed more accurately in the 1D model.
479!--       Update when opportunity arises!
480          IF ( humidity  .OR.  passive_scalar )  qs = 0.0
481
482!
483!--       inside buildings set velocities back to zero
484          IF ( topography /= 'flat' )  THEN
485             DO  i = nxl-1, nxr+1
486                DO  j = nys-1, nyn+1
487                   u(nzb:nzb_u_inner(j,i),j,i) = 0.0
488                   v(nzb:nzb_v_inner(j,i),j,i) = 0.0
489                ENDDO
490             ENDDO
491             IF ( conserve_volume_flow )  THEN
492                IF ( nxr == nx )  THEN
493                   DO  j = nys, nyn
494                      DO  k = nzb + 1, nzb_u_inner(j,nx)
495                         u_nzb_p1_for_vfc(j) = u_nzb_p1_for_vfc(j) + &
496                                               u1d(k) * dzu(k)
497                      ENDDO
498                   ENDDO
499                ENDIF
500                IF ( nyn == ny )  THEN
501                   DO  i = nxl, nxr
502                      DO  k = nzb + 1, nzb_v_inner(ny,i)
503                         v_nzb_p1_for_vfc(i) = v_nzb_p1_for_vfc(i) + &
504                                               v1d(k) * dzu(k)
505                      ENDDO
506                   ENDDO
507                ENDIF
508             ENDIF
509!
510!--          WARNING: The extra boundary conditions set after running the
511!--          -------  1D model impose an error on the divergence one layer
512!--                   below the topography; need to correct later
513!--          ATTENTION: Provisional correction for Piacsek & Williams
514!--          ---------  advection scheme: keep u and v zero one layer below
515!--                     the topography.
516             IF ( ibc_uv_b == 0 )  THEN
517!
518!--             Satisfying the Dirichlet condition with an extra layer below
519!--             the surface where the u and v component change their sign.
520                DO  i = nxl-1, nxr+1
521                   DO  j = nys-1, nyn+1
522                      IF ( nzb_u_inner(j,i) == 0 ) u(0,j,i) = -u(1,j,i)
523                      IF ( nzb_v_inner(j,i) == 0 ) v(0,j,i) = -v(1,j,i)
524                   ENDDO
525                ENDDO
526
527             ELSE
528!
529!--             Neumann condition
530                DO  i = nxl-1, nxr+1
531                   DO  j = nys-1, nyn+1
532                      IF ( nzb_u_inner(j,i) == 0 ) u(0,j,i) = u(1,j,i)
533                      IF ( nzb_v_inner(j,i) == 0 ) v(0,j,i) = v(1,j,i)
534                   ENDDO
535                ENDDO
536
537             ENDIF
538
539          ENDIF
540
541       ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 ) &
542       THEN
543!
544!--       Use constructed initial profiles (velocity constant with height,
545!--       temperature profile with constant gradient)
546          DO  i = nxl-1, nxr+1
547             DO  j = nys-1, nyn+1
548                pt(:,j,i) = pt_init
549                u(:,j,i)  = u_init
550                v(:,j,i)  = v_init
551             ENDDO
552          ENDDO
553
554!
555!--       Set initial horizontal velocities at the lowest computational grid
556!--       levels to zero in order to avoid too small time steps caused by the
557!--       diffusion limit in the initial phase of a run (at k=1, dz/2 occurs
558!--       in the limiting formula!). The original values are stored to be later
559!--       used for volume flow control.
560          DO  i = nxl-1, nxr+1
561             DO  j = nys-1, nyn+1
562                u(nzb:nzb_u_inner(j,i)+1,j,i) = 0.0
563                v(nzb:nzb_v_inner(j,i)+1,j,i) = 0.0
564             ENDDO
565          ENDDO
566          IF ( conserve_volume_flow )  THEN
567             IF ( nxr == nx )  THEN
568                DO  j = nys, nyn
569                   DO  k = nzb + 1, nzb_u_inner(j,nx) + 1
570                      u_nzb_p1_for_vfc(j) = u_nzb_p1_for_vfc(j) + &
571                                            u_init(k) * dzu(k)
572                   ENDDO
573                ENDDO
574             ENDIF
575             IF ( nyn == ny )  THEN
576                DO  i = nxl, nxr
577                   DO  k = nzb + 1, nzb_v_inner(ny,i) + 1
578                      v_nzb_p1_for_vfc(i) = v_nzb_p1_for_vfc(i) + &
579                                            v_init(k) * dzu(k)
580                   ENDDO
581                ENDDO
582             ENDIF
583          ENDIF
584
585          IF ( humidity  .OR.  passive_scalar )  THEN
586             DO  i = nxl-1, nxr+1
587                DO  j = nys-1, nyn+1
588                   q(:,j,i) = q_init
589                ENDDO
590             ENDDO
591          ENDIF
592
593          IF ( ocean )  THEN
594             DO  i = nxl-1, nxr+1
595                DO  j = nys-1, nyn+1
596                   sa(:,j,i) = sa_init
597                ENDDO
598             ENDDO
599          ENDIF
600         
601          IF ( constant_diffusion )  THEN
602             km   = km_constant
603             kh   = km / prandtl_number
604             e    = 0.0
605          ELSEIF ( e_init > 0.0 )  THEN
606             DO  k = nzb+1, nzt
607                km(k,:,:) = 0.1 * l_grid(k) * SQRT( e_init )
608             ENDDO
609             km(nzb,:,:)   = km(nzb+1,:,:)
610             km(nzt+1,:,:) = km(nzt,:,:)
611             kh   = km / prandtl_number
612             e    = e_init
613          ELSE
614             IF ( .NOT. ocean )  THEN
615                kh   = 0.01   ! there must exist an initial diffusion, because
616                km   = 0.01   ! otherwise no TKE would be produced by the
617                              ! production terms, as long as not yet
618                              ! e = (u*/cm)**2 at k=nzb+1
619             ELSE
620                kh   = 0.00001
621                km   = 0.00001
622             ENDIF
623             e    = 0.0
624          ENDIF
625          rif   = 0.0
626          ts    = 0.0
627          us    = 0.0
628          usws  = 0.0
629          uswst = top_momentumflux_u
630          vsws  = 0.0
631          vswst = top_momentumflux_v
632          IF ( humidity  .OR.  passive_scalar ) qs = 0.0
633
634!
635!--       Compute initial temperature field and other constants used in case
636!--       of a sloping surface
637          IF ( sloping_surface )  CALL init_slope
638
639       ELSEIF ( INDEX(initializing_actions, 'by_user') /= 0 ) &
640       THEN
641!
642!--       Initialization will completely be done by the user
643          CALL user_init_3d_model
644
645       ENDIF
646
647!
648!--    Apply channel flow boundary condition
649       IF ( TRIM( bc_uv_t ) == 'dirichlet_0' )  THEN
650
651          u(nzt+1,:,:) = 0.0
652          v(nzt+1,:,:) = 0.0
653
654!--       For the Dirichlet condition to be correctly applied at the top, set
655!--       ug and vg to zero there
656          ug(nzt+1)    = 0.0
657          vg(nzt+1)    = 0.0
658
659       ENDIF
660
661!
662!--    Calculate virtual potential temperature
663       IF ( humidity ) vpt = pt * ( 1.0 + 0.61 * q )
664
665!
666!--    Store initial profiles for output purposes etc.
667       hom(:,1,5,:) = SPREAD( u(:,nys,nxl), 2, statistic_regions+1 )
668       hom(:,1,6,:) = SPREAD( v(:,nys,nxl), 2, statistic_regions+1 )
669       IF ( ibc_uv_b == 0 )  THEN
670          hom(nzb,1,5,:) = -hom(nzb+1,1,5,:)  ! due to satisfying the Dirichlet
671          hom(nzb,1,6,:) = -hom(nzb+1,1,6,:)  ! condition with an extra layer
672              ! below the surface where the u and v component change their sign
673       ENDIF
674       hom(:,1,7,:)  = SPREAD( pt(:,nys,nxl), 2, statistic_regions+1 )
675       hom(:,1,23,:) = SPREAD( km(:,nys,nxl), 2, statistic_regions+1 )
676       hom(:,1,24,:) = SPREAD( kh(:,nys,nxl), 2, statistic_regions+1 )
677
678       IF ( ocean )  THEN
679!
680!--       Store initial salinity profile
681          hom(:,1,26,:)  = SPREAD( sa(:,nys,nxl), 2, statistic_regions+1 )
682       ENDIF
683
684       IF ( humidity )  THEN
685!
686!--       Store initial profile of total water content, virtual potential
687!--       temperature
688          hom(:,1,26,:) = SPREAD(   q(:,nys,nxl), 2, statistic_regions+1 )
689          hom(:,1,29,:) = SPREAD( vpt(:,nys,nxl), 2, statistic_regions+1 )
690          IF ( cloud_physics  .OR.  cloud_droplets ) THEN
691!
692!--          Store initial profile of specific humidity and potential
693!--          temperature
694             hom(:,1,27,:) = SPREAD(  q(:,nys,nxl), 2, statistic_regions+1 )
695             hom(:,1,28,:) = SPREAD( pt(:,nys,nxl), 2, statistic_regions+1 )
696          ENDIF
697       ENDIF
698
699       IF ( passive_scalar )  THEN
700!
701!--       Store initial scalar profile
702          hom(:,1,26,:) = SPREAD(  q(:,nys,nxl), 2, statistic_regions+1 )
703       ENDIF
704
705!
706!--    Initialize fluxes at bottom surface
707       IF ( use_surface_fluxes )  THEN
708
709          IF ( constant_heatflux )  THEN
710!
711!--          Heat flux is prescribed
712             IF ( random_heatflux )  THEN
713                CALL disturb_heatflux
714             ELSE
715                shf = surface_heatflux
716!
717!--             Over topography surface_heatflux is replaced by wall_heatflux(0)
718                IF ( TRIM( topography ) /= 'flat' )  THEN
719                   DO  i = nxl-1, nxr+1
720                      DO  j = nys-1, nyn+1
721                         IF ( nzb_s_inner(j,i) /= 0 )  THEN
722                            shf(j,i) = wall_heatflux(0)
723                         ENDIF
724                      ENDDO
725                   ENDDO
726                ENDIF
727             ENDIF
728             IF ( ASSOCIATED( shf_m ) )  shf_m = shf
729          ENDIF
730
731!
732!--       Determine the near-surface water flux
733          IF ( humidity  .OR.  passive_scalar )  THEN
734             IF ( constant_waterflux )  THEN
735                qsws   = surface_waterflux
736                IF ( ASSOCIATED( qsws_m ) )  qsws_m = qsws
737             ENDIF
738          ENDIF
739
740       ENDIF
741
742!
743!--    Initialize fluxes at top surface
744!--    Currently, only the heatflux and salinity flux can be prescribed.
745!--    The latent flux is zero in this case!
746       IF ( use_top_fluxes )  THEN
747
748          IF ( constant_top_heatflux )  THEN
749!
750!--          Heat flux is prescribed
751             tswst = top_heatflux
752             IF ( ASSOCIATED( tswst_m ) )  tswst_m = tswst
753
754             IF ( humidity  .OR.  passive_scalar )  THEN
755                qswst = 0.0
756                IF ( ASSOCIATED( qswst_m ) )  qswst_m = qswst
757             ENDIF
758
759             IF ( ocean )  THEN
760                saswsb = bottom_salinityflux
761                saswst = top_salinityflux
762             ENDIF
763          ENDIF
764
765!
766!--       Initialization in case of a coupled model run
767          IF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
768             tswst = 0.0
769             IF ( ASSOCIATED( tswst_m ) )  tswst_m = tswst
770          ENDIF
771
772       ENDIF
773
774!
775!--    Initialize Prandtl layer quantities
776       IF ( prandtl_layer )  THEN
777
778          z0 = roughness_length
779
780          IF ( .NOT. constant_heatflux )  THEN 
781!
782!--          Surface temperature is prescribed. Here the heat flux cannot be
783!--          simply estimated, because therefore rif, u* and theta* would have
784!--          to be computed by iteration. This is why the heat flux is assumed
785!--          to be zero before the first time step. It approaches its correct
786!--          value in the course of the first few time steps.
787             shf   = 0.0
788             IF ( ASSOCIATED( shf_m ) )  shf_m = 0.0
789          ENDIF
790
791          IF ( humidity  .OR.  passive_scalar )  THEN
792             IF ( .NOT. constant_waterflux )  THEN
793                qsws   = 0.0
794                IF ( ASSOCIATED( qsws_m ) )   qsws_m = 0.0
795             ENDIF
796          ENDIF
797
798       ENDIF
799
800!
801!--    Calculate the initial volume flow at the right and north boundary
802       IF ( conserve_volume_flow )  THEN
803
804          volume_flow_initial_l = 0.0
805          volume_flow_area_l    = 0.0
806 
807          IF ( nxr == nx )  THEN
808             DO  j = nys, nyn
809                DO  k = nzb_2d(j,nx) + 1, nzt
810                   volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
811                                              u(k,j,nx) * dzu(k)
812                   volume_flow_area_l(1)    = volume_flow_area_l(1) + dzu(k)
813                ENDDO
814!
815!--             Correction if velocity at nzb+1 has been set zero further above
816                volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
817                                           u_nzb_p1_for_vfc(j)
818             ENDDO
819          ENDIF
820
821          IF ( nyn == ny )  THEN
822             DO  i = nxl, nxr
823                DO  k = nzb_2d(ny,i) + 1, nzt 
824                   volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
825                                              v(k,ny,i) * dzu(k)
826                   volume_flow_area_l(2)    = volume_flow_area_l(2) + dzu(k)
827                ENDDO
828!
829!--             Correction if velocity at nzb+1 has been set zero further above
830                volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
831                                           v_nzb_p1_for_vfc(i)
832             ENDDO
833          ENDIF
834
835#if defined( __parallel )
836          CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&
837                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
838          CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1),      &
839                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
840#else
841          volume_flow_initial = volume_flow_initial_l
842          volume_flow_area    = volume_flow_area_l
843#endif
844!
845!--       In case of 'bulk_velocity' mode, volume_flow_initial is overridden
846!--       and calculated from u|v_bulk instead.
847          IF ( TRIM( conserve_volume_flow_mode ) == 'bulk_velocity' )  THEN
848             volume_flow_initial(1) = u_bulk * volume_flow_area(1)
849             volume_flow_initial(2) = v_bulk * volume_flow_area(2)
850          ENDIF
851
852       ENDIF
853
854!
855!--    For the moment, perturbation pressure and vertical velocity are zero
856       p = 0.0; w = 0.0
857
858!
859!--    Initialize array sums (must be defined in first call of pres)
860       sums = 0.0
861
862!
863!--    Treating cloud physics, liquid water content and precipitation amount
864!--    are zero at beginning of the simulation
865       IF ( cloud_physics )  THEN
866          ql = 0.0
867          IF ( precipitation )  precipitation_amount = 0.0
868       ENDIF
869
870!
871!--    Impose vortex with vertical axis on the initial velocity profile
872       IF ( INDEX( initializing_actions, 'initialize_vortex' ) /= 0 )  THEN
873          CALL init_rankine
874       ENDIF
875
876!
877!--    Impose temperature anomaly (advection test only)
878       IF ( INDEX( initializing_actions, 'initialize_ptanom' ) /= 0 )  THEN
879          CALL init_pt_anomaly
880       ENDIF
881
882!
883!--    If required, change the surface temperature at the start of the 3D run
884       IF ( pt_surface_initial_change /= 0.0 )  THEN
885          pt(nzb,:,:) = pt(nzb,:,:) + pt_surface_initial_change
886       ENDIF
887
888!
889!--    If required, change the surface humidity/scalar at the start of the 3D
890!--    run
891       IF ( ( humidity .OR. passive_scalar ) .AND. &
892            q_surface_initial_change /= 0.0 )  THEN
893          q(nzb,:,:) = q(nzb,:,:) + q_surface_initial_change
894       ENDIF
895
896!
897!--    Initialize the random number generator (from numerical recipes)
898       CALL random_function_ini
899
900!
901!--    Impose random perturbation on the horizontal velocity field and then
902!--    remove the divergences from the velocity field
903       IF ( create_disturbances )  THEN
904          CALL disturb_field( nzb_u_inner, tend, u )
905          CALL disturb_field( nzb_v_inner, tend, v )
906          n_sor = nsor_ini
907          CALL pres
908          n_sor = nsor
909       ENDIF
910
911!
912!--    Once again set the perturbation pressure explicitly to zero in order to
913!--    assure that it does not generate any divergences in the first time step.
914!--    At t=0 the velocity field is free of divergence (as constructed above).
915!--    Divergences being created during a time step are not yet known and thus
916!--    cannot be corrected during the time step yet.
917       p = 0.0
918
919!
920!--    Initialize old and new time levels.
921       IF ( timestep_scheme(1:5) /= 'runge' )  THEN
922          e_m = e; pt_m = pt; u_m = u; v_m = v; w_m = w; kh_m = kh; km_m = km
923       ELSE
924          te_m = 0.0; tpt_m = 0.0; tu_m = 0.0; tv_m = 0.0; tw_m = 0.0
925       ENDIF
926       e_p = e; pt_p = pt; u_p = u; v_p = v; w_p = w
927
928       IF ( humidity  .OR.  passive_scalar )  THEN
929          IF ( ASSOCIATED( q_m ) )  q_m = q
930          IF ( timestep_scheme(1:5) == 'runge' )  tq_m = 0.0
931          q_p = q
932          IF ( humidity  .AND.  ASSOCIATED( vpt_m ) )  vpt_m = vpt
933       ENDIF
934
935       IF ( ocean )  THEN
936          tsa_m = 0.0
937          sa_p  = sa
938       ENDIF
939
940
941    ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data'  .OR.    &
942             TRIM( initializing_actions ) == 'cyclic_fill' )  &
943    THEN
944!
945!--    When reading data for cyclic fill of 3D prerun data, first read
946!--    some of the global variables from restart file
947       IF ( TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
948
949          WRITE (9,*) 'before read_parts_of_var_list'
950          CALL local_flush( 9 )
951          CALL read_parts_of_var_list
952          WRITE (9,*) 'after read_parts_of_var_list'
953          CALL local_flush( 9 )
954          CALL close_file( 13 )
955
956!
957!--       Initialization of the turbulence recycling method
958          IF ( turbulent_inflow )  THEN
959!
960!--          Store temporally and horizontally averaged vertical profiles to be
961!--          used as mean inflow profiles
962             ALLOCATE( mean_inflow_profiles(nzb:nzt+1,5) )
963
964             mean_inflow_profiles(:,1) = hom_sum(:,1,0)    ! u
965             mean_inflow_profiles(:,2) = hom_sum(:,2,0)    ! v
966             mean_inflow_profiles(:,4) = hom_sum(:,4,0)    ! pt
967             mean_inflow_profiles(:,5) = hom_sum(:,8,0)    ! e
968
969!
970!--          Use these mean profiles for the inflow (provided that Dirichlet
971!--          conditions are used)
972             IF ( inflow_l )  THEN
973                DO  j = nys-1, nyn+1
974                   DO  k = nzb, nzt+1
975                      u(k,j,-1)  = mean_inflow_profiles(k,1)
976                      v(k,j,-1)  = mean_inflow_profiles(k,2)
977                      w(k,j,-1)  = 0.0
978                      pt(k,j,-1) = mean_inflow_profiles(k,4)
979                      e(k,j,-1)  = mean_inflow_profiles(k,5)
980                   ENDDO
981                ENDDO
982             ENDIF
983
984!
985!--          Calculate the damping factors to be used at the inflow. For a
986!--          turbulent inflow the turbulent fluctuations have to be limited
987!--          vertically because otherwise the turbulent inflow layer will grow
988!--          in time.
989             IF ( inflow_damping_height == 9999999.9 )  THEN
990!
991!--             Default: use the inversion height calculated by the prerun; if
992!--             this is zero, inflow_damping_height must be explicitly
993!--             specified.
994                IF ( hom_sum(nzb+6,pr_palm,0) /= 0.0 )  THEN
995                   inflow_damping_height = hom_sum(nzb+6,pr_palm,0)
996                ELSE
997                   WRITE( message_string, * ) 'inflow_damping_height must be ',&
998                        'explicitly specified because&the inversion height ', &
999                        'calculated by the prerun is zero.'
1000                   CALL message( 'init_3d_model', 'PA0318', 1, 2, 0, 6, 0 )
1001                ENDIF
1002
1003             ENDIF
1004
1005             IF ( inflow_damping_width == 9999999.9 )  THEN
1006!
1007!--             Default for the transition range: one tenth of the undamped
1008!--             layer
1009                inflow_damping_width = 0.1 * inflow_damping_height
1010
1011             ENDIF
1012
1013             ALLOCATE( inflow_damping_factor(nzb:nzt+1) )
1014
1015             DO  k = nzb, nzt+1
1016
1017                IF ( zu(k) <= inflow_damping_height )  THEN
1018                   inflow_damping_factor(k) = 1.0
1019                ELSEIF ( zu(k) <= inflow_damping_height +  &
1020                                  inflow_damping_width )  THEN
1021                   inflow_damping_factor(k) = 1.0 -                            &
1022                                           ( zu(k) - inflow_damping_height ) / &
1023                                           inflow_damping_width
1024                ELSE
1025                   inflow_damping_factor(k) = 0.0
1026                ENDIF
1027
1028             ENDDO
1029          ENDIF
1030
1031       ENDIF
1032
1033!
1034!--    Read binary data from restart file
1035          WRITE (9,*) 'before read_3d_binary'
1036          CALL local_flush( 9 )
1037       CALL read_3d_binary
1038          WRITE (9,*) 'after read_3d_binary'
1039          CALL local_flush( 9 )
1040
1041!
1042!--    Inside buildings set velocities and TKE back to zero
1043       IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND.  &
1044            topography /= 'flat' )  THEN
1045!
1046!--       Correction of initial volume flow
1047          IF ( conserve_volume_flow )  THEN
1048             IF ( nxr == nx )  THEN
1049                DO  j = nys, nyn
1050                   DO  k = nzb + 1, nzb_u_inner(j,nx)
1051                      u_nzb_p1_for_vfc(j) = u_nzb_p1_for_vfc(j) + &
1052                                            u(k,j,nx) * dzu(k)
1053                   ENDDO
1054                ENDDO
1055             ENDIF
1056             IF ( nyn == ny )  THEN
1057                DO  i = nxl, nxr
1058                   DO  k = nzb + 1, nzb_v_inner(ny,i)
1059                      v_nzb_p1_for_vfc(i) = v_nzb_p1_for_vfc(i) + &
1060                                            v(k,ny,i) * dzu(k)
1061                   ENDDO
1062                ENDDO
1063             ENDIF
1064          ENDIF
1065
1066!
1067!--       Inside buildings set velocities and TKE back to zero.
1068!--       Other scalars (pt, q, s, km, kh, p, sa, ...) are ignored at present,
1069!--       maybe revise later.
1070          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1071             DO  i = nxl-1, nxr+1
1072                DO  j = nys-1, nyn+1
1073                   u  (nzb:nzb_u_inner(j,i),j,i) = 0.0
1074                   v  (nzb:nzb_v_inner(j,i),j,i) = 0.0
1075                   w  (nzb:nzb_w_inner(j,i),j,i) = 0.0
1076                   e  (nzb:nzb_w_inner(j,i),j,i) = 0.0
1077                   u_m(nzb:nzb_u_inner(j,i),j,i) = 0.0
1078                   v_m(nzb:nzb_v_inner(j,i),j,i) = 0.0
1079                   w_m(nzb:nzb_w_inner(j,i),j,i) = 0.0
1080                   e_m(nzb:nzb_w_inner(j,i),j,i) = 0.0
1081                   tu_m(nzb:nzb_u_inner(j,i),j,i) = 0.0
1082                   tv_m(nzb:nzb_v_inner(j,i),j,i) = 0.0
1083                   tw_m(nzb:nzb_w_inner(j,i),j,i) = 0.0
1084                   te_m(nzb:nzb_w_inner(j,i),j,i) = 0.0
1085                   tpt_m(nzb:nzb_w_inner(j,i),j,i) = 0.0
1086                ENDDO
1087             ENDDO
1088          ELSE
1089             DO  i = nxl-1, nxr+1
1090                DO  j = nys-1, nyn+1
1091                   u  (nzb:nzb_u_inner(j,i),j,i) = 0.0
1092                   v  (nzb:nzb_v_inner(j,i),j,i) = 0.0
1093                   w  (nzb:nzb_w_inner(j,i),j,i) = 0.0
1094                   e  (nzb:nzb_w_inner(j,i),j,i) = 0.0
1095                   u_m(nzb:nzb_u_inner(j,i),j,i) = 0.0
1096                   v_m(nzb:nzb_v_inner(j,i),j,i) = 0.0
1097                   w_m(nzb:nzb_w_inner(j,i),j,i) = 0.0
1098                   e_m(nzb:nzb_w_inner(j,i),j,i) = 0.0
1099                   u_p(nzb:nzb_u_inner(j,i),j,i) = 0.0
1100                   v_p(nzb:nzb_v_inner(j,i),j,i) = 0.0
1101                   w_p(nzb:nzb_w_inner(j,i),j,i) = 0.0
1102                   e_p(nzb:nzb_w_inner(j,i),j,i) = 0.0
1103                ENDDO
1104             ENDDO
1105          ENDIF
1106
1107       ENDIF
1108
1109!
1110!--    Calculate the initial volume flow at the right and north boundary
1111       IF ( conserve_volume_flow  .AND.  &
1112            TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
1113
1114          volume_flow_initial_l = 0.0
1115          volume_flow_area_l    = 0.0
1116 
1117          IF ( nxr == nx )  THEN
1118             DO  j = nys, nyn
1119                DO  k = nzb_2d(j,nx) + 1, nzt
1120                   volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
1121                                              u(k,j,nx) * dzu(k)
1122                   volume_flow_area_l(1)    = volume_flow_area_l(1) + dzu(k)
1123                ENDDO
1124!
1125!--             Correction if velocity inside buildings has been set to zero
1126!--             further above
1127                volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
1128                                           u_nzb_p1_for_vfc(j)
1129             ENDDO
1130          ENDIF
1131
1132          IF ( nyn == ny )  THEN
1133             DO  i = nxl, nxr
1134                DO  k = nzb_2d(ny,i) + 1, nzt 
1135                   volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
1136                                              v(k,ny,i) * dzu(k)
1137                   volume_flow_area_l(2)    = volume_flow_area_l(2) + dzu(k)
1138                ENDDO
1139!
1140!--             Correction if velocity inside buildings has been set to zero
1141!--             further above
1142                volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
1143                                           v_nzb_p1_for_vfc(i)
1144             ENDDO
1145          ENDIF
1146
1147#if defined( __parallel )
1148          CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&
1149                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
1150          CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1),      &
1151                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
1152#else
1153          volume_flow_initial = volume_flow_initial_l
1154          volume_flow_area    = volume_flow_area_l
1155#endif 
1156       ENDIF
1157
1158
1159!
1160!--    Calculate initial temperature field and other constants used in case
1161!--    of a sloping surface
1162       IF ( sloping_surface )  CALL init_slope
1163
1164!
1165!--    Initialize new time levels (only done in order to set boundary values
1166!--    including ghost points)
1167       e_p = e; pt_p = pt; u_p = u; v_p = v; w_p = w
1168       IF ( humidity  .OR.  passive_scalar )  q_p = q
1169       IF ( ocean )  sa_p = sa
1170
1171!
1172!--    Allthough tendency arrays are set in prognostic_equations, they have
1173!--    have to be predefined here because they are used (but multiplied with 0)
1174!--    there before they are set.
1175       IF ( timestep_scheme(1:5) == 'runge' )  THEN
1176          te_m = 0.0; tpt_m = 0.0; tu_m = 0.0; tv_m = 0.0; tw_m = 0.0
1177          IF ( humidity  .OR.  passive_scalar )  tq_m = 0.0
1178          IF ( ocean )  tsa_m = 0.0
1179       ENDIF
1180
1181    ELSE
1182!
1183!--    Actually this part of the programm should not be reached
1184       message_string = 'unknown initializing problem'
1185       CALL message( 'init_3d_model', 'PA0193', 1, 2, 0, 6, 0 )
1186    ENDIF
1187
1188
1189    IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1190!
1191!--    Initialize old timelevels needed for radiation boundary conditions
1192       IF ( outflow_l )  THEN
1193          u_m_l(:,:,:) = u(:,:,1:2)
1194          v_m_l(:,:,:) = v(:,:,0:1)
1195          w_m_l(:,:,:) = w(:,:,0:1)
1196       ENDIF
1197       IF ( outflow_r )  THEN
1198          u_m_r(:,:,:) = u(:,:,nx-1:nx)
1199          v_m_r(:,:,:) = v(:,:,nx-1:nx)
1200          w_m_r(:,:,:) = w(:,:,nx-1:nx)
1201       ENDIF
1202       IF ( outflow_s )  THEN
1203          u_m_s(:,:,:) = u(:,0:1,:)
1204          v_m_s(:,:,:) = v(:,1:2,:)
1205          w_m_s(:,:,:) = w(:,0:1,:)
1206       ENDIF
1207       IF ( outflow_n )  THEN
1208          u_m_n(:,:,:) = u(:,ny-1:ny,:)
1209          v_m_n(:,:,:) = v(:,ny-1:ny,:)
1210          w_m_n(:,:,:) = w(:,ny-1:ny,:)
1211       ENDIF
1212
1213    ENDIF
1214
1215!
1216!-- Initialization of the leaf area density
1217    IF ( plant_canopy ) THEN
1218 
1219       SELECT CASE ( TRIM( canopy_mode ) )
1220
1221          CASE( 'block' )
1222
1223             DO  i = nxl-1, nxr+1
1224                DO  j = nys-1, nyn+1
1225                   lad_s(:,j,i) = lad(:)
1226                   cdc(:,j,i)   = drag_coefficient
1227                   IF ( passive_scalar ) THEN
1228                      sls(:,j,i) = leaf_surface_concentration
1229                      sec(:,j,i) = scalar_exchange_coefficient
1230                   ENDIF
1231                ENDDO
1232             ENDDO
1233
1234          CASE DEFAULT
1235
1236!
1237!--          The DEFAULT case is reached either if the parameter
1238!--          canopy mode contains a wrong character string or if the
1239!--          user has coded a special case in the user interface.
1240!--          There, the subroutine user_init_plant_canopy checks
1241!--          which of these two conditions applies.
1242             CALL user_init_plant_canopy
1243 
1244          END SELECT
1245
1246       CALL exchange_horiz( lad_s )
1247       CALL exchange_horiz( cdc )
1248
1249       IF ( passive_scalar ) THEN
1250          CALL exchange_horiz( sls )
1251          CALL exchange_horiz( sec )
1252       ENDIF
1253
1254!
1255!--    Sharp boundaries of the plant canopy in horizontal directions
1256!--    In vertical direction the interpolation is retained, as the leaf
1257!--    area density is initialised by prescribing a vertical profile
1258!--    consisting of piecewise linear segments. The upper boundary
1259!--    of the plant canopy is now defined by lad_w(pch_index,:,:) = 0.0.
1260
1261       DO  i = nxl, nxr
1262          DO  j = nys, nyn
1263             DO  k = nzb, nzt+1 
1264                IF ( lad_s(k,j,i) > 0.0 ) THEN
1265                   lad_u(k,j,i)   = lad_s(k,j,i) 
1266                   lad_u(k,j,i+1) = lad_s(k,j,i)
1267                   lad_v(k,j,i)   = lad_s(k,j,i)
1268                   lad_v(k,j+1,i) = lad_s(k,j,i)
1269                ENDIF
1270             ENDDO
1271             DO  k = nzb, nzt
1272                lad_w(k,j,i) = 0.5 * ( lad_s(k+1,j,i) + lad_s(k,j,i) )
1273             ENDDO
1274          ENDDO
1275       ENDDO
1276
1277       lad_w(pch_index,:,:) = 0.0
1278       lad_w(nzt+1,:,:)     = lad_w(nzt,:,:)
1279
1280       CALL exchange_horiz( lad_u )
1281       CALL exchange_horiz( lad_v )
1282       CALL exchange_horiz( lad_w )
1283
1284!
1285!--    Initialisation of the canopy heat source distribution
1286       IF ( cthf /= 0.0 ) THEN
1287!
1288!--       Piecewise evaluation of the leaf area index by
1289!--       integration of the leaf area density
1290          lai(:,:,:) = 0.0
1291          DO  i = nxl-1, nxr+1
1292             DO  j = nys-1, nyn+1
1293                DO  k = pch_index-1, 0, -1
1294                   lai(k,j,i) = lai(k+1,j,i) +                   &
1295                                ( 0.5 * ( lad_w(k+1,j,i) +       &
1296                                          lad_s(k+1,j,i) ) *     &
1297                                  ( zw(k+1) - zu(k+1) ) )  +     &
1298                                ( 0.5 * ( lad_w(k,j,i)   +       &
1299                                          lad_s(k+1,j,i) ) *     &
1300                                  ( zu(k+1) - zw(k) ) )
1301                ENDDO
1302             ENDDO
1303          ENDDO
1304
1305!
1306!--       Evaluation of the upward kinematic vertical heat flux within the
1307!--       canopy
1308          DO  i = nxl-1, nxr+1
1309             DO  j = nys-1, nyn+1
1310                DO  k = 0, pch_index
1311                   canopy_heat_flux(k,j,i) = cthf *                    &
1312                                             exp( -0.6 * lai(k,j,i) )
1313                ENDDO
1314             ENDDO
1315          ENDDO
1316
1317!
1318!--       The near surface heat flux is derived from the heat flux
1319!--       distribution within the canopy
1320          shf(:,:) = canopy_heat_flux(0,:,:)
1321
1322          IF ( ASSOCIATED( shf_m ) ) shf_m = shf
1323
1324       ENDIF
1325
1326    ENDIF
1327
1328!
1329!-- If required, initialize dvrp-software
1330    IF ( dt_dvrp /= 9999999.9 )  CALL init_dvrp
1331
1332    IF ( ocean )  THEN
1333!
1334!--    Initialize quantities needed for the ocean model
1335       CALL init_ocean
1336
1337    ELSE
1338!
1339!--    Initialize quantities for handling cloud physics
1340!--    This routine must be called before init_particles, because
1341!--    otherwise, array pt_d_t, needed in data_output_dvrp (called by
1342!--    init_particles) is not defined.
1343       CALL init_cloud_physics
1344    ENDIF
1345
1346!
1347!-- If required, initialize particles
1348    IF ( particle_advection )  CALL init_particles
1349
1350!
1351!-- Initialize quantities for special advections schemes
1352    CALL init_advec
1353
1354!
1355!-- Initialize Rayleigh damping factors
1356    rdf = 0.0
1357    IF ( rayleigh_damping_factor /= 0.0 )  THEN
1358       IF ( .NOT. ocean )  THEN
1359          DO  k = nzb+1, nzt
1360             IF ( zu(k) >= rayleigh_damping_height )  THEN
1361                rdf(k) = rayleigh_damping_factor * &
1362                      ( SIN( pi * 0.5 * ( zu(k) - rayleigh_damping_height )    &
1363                                      / ( zu(nzt) - rayleigh_damping_height ) )&
1364                      )**2
1365             ENDIF
1366          ENDDO
1367       ELSE
1368          DO  k = nzt, nzb+1, -1
1369             IF ( zu(k) <= rayleigh_damping_height )  THEN
1370                rdf(k) = rayleigh_damping_factor * &
1371                      ( SIN( pi * 0.5 * ( rayleigh_damping_height - zu(k) )    &
1372                                      / ( rayleigh_damping_height - zu(nzb+1)))&
1373                      )**2
1374             ENDIF
1375          ENDDO
1376       ENDIF
1377    ENDIF
1378
1379!
1380!-- Initialize the starting level and the vertical smoothing factor used for
1381!-- the external pressure gradient
1382    dp_smooth_factor = 1.0
1383    IF ( dp_external )  THEN
1384!
1385!--    Set the starting level dp_level_ind_b only if it has not been set before
1386!--    (e.g. in init_grid).
1387       IF ( dp_level_ind_b == 0 )  THEN
1388          ind_array = MINLOC( ABS( dp_level_b - zu ) )
1389          dp_level_ind_b = ind_array(1) - 1 + nzb
1390                                        ! MINLOC uses lower array bound 1
1391       ENDIF
1392       IF ( dp_smooth )  THEN
1393          dp_smooth_factor(:dp_level_ind_b) = 0.0
1394          DO  k = dp_level_ind_b+1, nzt
1395             dp_smooth_factor(k) = 0.5 * ( 1.0 + SIN( pi * &
1396                  ( REAL( k - dp_level_ind_b ) /  &
1397                    REAL( nzt - dp_level_ind_b ) - 0.5 ) ) )
1398          ENDDO
1399       ENDIF
1400    ENDIF
1401
1402!
1403!-- Initialize diffusivities used within the outflow damping layer in case of
1404!-- non-cyclic lateral boundaries. A linear increase is assumed over the first
1405!-- half of the width of the damping layer
1406    IF ( bc_lr == 'dirichlet/radiation' )  THEN
1407
1408       DO  i = nxl-1, nxr+1
1409          IF ( i >= nx - outflow_damping_width )  THEN
1410             km_damp_x(i) = km_damp_max * MIN( 1.0,                    &
1411                            ( i - ( nx - outflow_damping_width ) ) /   &
1412                            REAL( outflow_damping_width/2 )            &
1413                                             )
1414          ELSE
1415             km_damp_x(i) = 0.0
1416          ENDIF
1417       ENDDO
1418
1419    ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
1420
1421       DO  i = nxl-1, nxr+1
1422          IF ( i <= outflow_damping_width )  THEN
1423             km_damp_x(i) = km_damp_max * MIN( 1.0,                    &
1424                            ( outflow_damping_width - i ) /            &
1425                            REAL( outflow_damping_width/2 )            &
1426                                             )
1427          ELSE
1428             km_damp_x(i) = 0.0
1429          ENDIF
1430       ENDDO
1431
1432    ENDIF
1433
1434    IF ( bc_ns == 'dirichlet/radiation' )  THEN
1435
1436       DO  j = nys-1, nyn+1
1437          IF ( j >= ny - outflow_damping_width )  THEN
1438             km_damp_y(j) = km_damp_max * MIN( 1.0,                    &
1439                            ( j - ( ny - outflow_damping_width ) ) /   &
1440                            REAL( outflow_damping_width/2 )            &
1441                                             )
1442          ELSE
1443             km_damp_y(j) = 0.0
1444          ENDIF
1445       ENDDO
1446
1447    ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
1448
1449       DO  j = nys-1, nyn+1
1450          IF ( j <= outflow_damping_width )  THEN
1451             km_damp_y(j) = km_damp_max * MIN( 1.0,                    &
1452                            ( outflow_damping_width - j ) /            &
1453                            REAL( outflow_damping_width/2 )            &
1454                                             )
1455          ELSE
1456             km_damp_y(j) = 0.0
1457          ENDIF
1458       ENDDO
1459
1460    ENDIF
1461
1462!
1463!-- Initialize local summation arrays for UP flow_statistics. This is necessary
1464!-- because they may not yet have been initialized when they are called from
1465!-- flow_statistics (or - depending on the chosen model run - are never
1466!-- initialized)
1467    sums_divnew_l      = 0.0
1468    sums_divold_l      = 0.0
1469    sums_l_l           = 0.0
1470    sums_up_fraction_l = 0.0
1471    sums_wsts_bc_l     = 0.0
1472
1473!
1474!-- Pre-set masks for regional statistics. Default is the total model domain.
1475    rmask = 1.0
1476
1477!
1478!-- User-defined initializing actions. Check afterwards, if maximum number
1479!-- of allowed timeseries is not exceeded
1480    CALL user_init
1481
1482    IF ( dots_num > dots_max )  THEN
1483       WRITE( message_string, * ) 'number of time series quantities exceeds', &
1484                                  ' its maximum of dots_max = ', dots_max,    &
1485                                  ' &Please increase dots_max in modules.f90.'
1486       CALL message( 'init_3d_model', 'PA0194', 1, 2, 0, 6, 0 )   
1487    ENDIF
1488
1489!
1490!-- Input binary data file is not needed anymore. This line must be placed
1491!-- after call of user_init!
1492    CALL close_file( 13 )
1493
1494!
1495!-- Compute total sum of active mask grid points
1496!-- ngp_2dh: number of grid points of a horizontal cross section through the
1497!--          total domain
1498!-- ngp_3d:  number of grid points of the total domain
1499    ngp_2dh_outer_l   = 0
1500    ngp_2dh_outer     = 0
1501    ngp_2dh_s_inner_l = 0
1502    ngp_2dh_s_inner   = 0
1503    ngp_2dh_l         = 0
1504    ngp_2dh           = 0
1505    ngp_3d_inner_l    = 0
1506    ngp_3d_inner      = 0
1507    ngp_3d            = 0
1508    ngp_sums          = ( nz + 2 ) * ( pr_palm + max_pr_user )
1509
1510    DO  sr = 0, statistic_regions
1511       DO  i = nxl, nxr
1512          DO  j = nys, nyn
1513             IF ( rmask(j,i,sr) == 1.0 )  THEN
1514!
1515!--             All xy-grid points
1516                ngp_2dh_l(sr) = ngp_2dh_l(sr) + 1
1517!
1518!--             xy-grid points above topography
1519                DO  k = nzb_s_outer(j,i), nz + 1
1520                   ngp_2dh_outer_l(k,sr) = ngp_2dh_outer_l(k,sr) + 1
1521                ENDDO
1522                DO  k = nzb_s_inner(j,i), nz + 1
1523                   ngp_2dh_s_inner_l(k,sr) = ngp_2dh_s_inner_l(k,sr) + 1
1524                ENDDO
1525!
1526!--             All grid points of the total domain above topography
1527                ngp_3d_inner_l(sr) = ngp_3d_inner_l(sr) + &
1528                                     ( nz - nzb_s_inner(j,i) + 2 )
1529             ENDIF
1530          ENDDO
1531       ENDDO
1532    ENDDO
1533
1534    sr = statistic_regions + 1
1535#if defined( __parallel )
1536    CALL MPI_ALLREDUCE( ngp_2dh_l(0), ngp_2dh(0), sr, MPI_INTEGER, MPI_SUM,  &
1537                        comm2d, ierr )
1538    CALL MPI_ALLREDUCE( ngp_2dh_outer_l(0,0), ngp_2dh_outer(0,0), (nz+2)*sr, &
1539                        MPI_INTEGER, MPI_SUM, comm2d, ierr )
1540    CALL MPI_ALLREDUCE( ngp_2dh_s_inner_l(0,0), ngp_2dh_s_inner(0,0),        &
1541                        (nz+2)*sr, MPI_INTEGER, MPI_SUM, comm2d, ierr )
1542    CALL MPI_ALLREDUCE( ngp_3d_inner_l(0), ngp_3d_inner(0), sr, MPI_INTEGER, &
1543                        MPI_SUM, comm2d, ierr )
1544#else
1545    ngp_2dh         = ngp_2dh_l
1546    ngp_2dh_outer   = ngp_2dh_outer_l
1547    ngp_2dh_s_inner = ngp_2dh_s_inner_l
1548    ngp_3d_inner    = ngp_3d_inner_l
1549#endif
1550
1551    ngp_3d = ngp_2dh * ( nz + 2 )
1552
1553!
1554!-- Set a lower limit of 1 in order to avoid zero divisions in flow_statistics,
1555!-- buoyancy, etc. A zero value will occur for cases where all grid points of
1556!-- the respective subdomain lie below the surface topography
1557    ngp_2dh_outer   = MAX( 1, ngp_2dh_outer(:,:)   ) 
1558    ngp_3d_inner    = MAX( 1, ngp_3d_inner(:)      )
1559    ngp_2dh_s_inner = MAX( 1, ngp_2dh_s_inner(:,:) ) 
1560
1561    DEALLOCATE( ngp_2dh_l, ngp_2dh_outer_l, ngp_3d_inner_l )
1562
1563
1564 END SUBROUTINE init_3d_model
Note: See TracBrowser for help on using the repository browser.