source: palm/trunk/SOURCE/init_3d_model.f90 @ 2007

Last change on this file since 2007 was 2007, checked in by kanani, 5 years ago

changes in the course of urban surface model implementation

  • Property svn:keywords set to Id
  • Property svn:mergeinfo set to False
    /palm/branches/forwind/SOURCE/init_3d_model.f901564-1913
File size: 74.8 KB
Line 
1!> @file init_3d_model.f90
2!------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2016 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22! Added support for urban surface model,
23! adjusted location_message in case of plant_canopy
24!
25! Former revisions:
26! -----------------
27! $Id: init_3d_model.f90 2007 2016-08-24 15:47:17Z kanani $
28!
29! 2000 2016-08-20 18:09:15Z knoop
30! Forced header and separation lines into 80 columns
31!
32! 1992 2016-08-12 15:14:59Z suehring
33! Initializaton of scalarflux at model top
34! Bugfixes in initialization of surface and top salinity flux, top scalar and
35! humidity fluxes
36!
37! 1960 2016-07-12 16:34:24Z suehring
38! Separate humidity and passive scalar
39! Increase dimension for mean_inflow_profiles
40! Remove inadvertent write-statement
41! Bugfix, large-scale forcing is still not implemented for passive scalars
42!
43! 1957 2016-07-07 10:43:48Z suehring
44! flight module added
45!
46! 1920 2016-05-30 10:50:15Z suehring
47! Initialize us with very small number to avoid segmentation fault during
48! calculation of Obukhov length
49!
50! 1918 2016-05-27 14:35:57Z raasch
51! intermediate_timestep_count is set 0 instead 1 for first call of pres,
52! bugfix: initialization of local sum arrays are moved to the beginning of the
53!         routine because otherwise results from pres are overwritten
54!
55! 1914 2016-05-26 14:44:07Z witha
56! Added initialization of the wind turbine model
57!
58! 1878 2016-04-19 12:30:36Z hellstea
59! The zeroth element of weight_pres removed as unnecessary
60!
61! 1849 2016-04-08 11:33:18Z hoffmann
62! Adapted for modularization of microphysics.
63! precipitation_amount, precipitation_rate, prr moved to arrays_3d.
64! Initialization of nc_1d, nr_1d, pt_1d, qc_1d, qr_1d, q_1d moved to
65! microphysics_init.
66!
67! 1845 2016-04-08 08:29:13Z raasch
68! nzb_2d replaced by nzb_u|v_inner
69!
70! 1833 2016-04-07 14:23:03Z raasch
71! initialization of spectra quantities moved to spectra_mod
72!
73! 1831 2016-04-07 13:15:51Z hoffmann
74! turbulence renamed collision_turbulence
75!
76! 1826 2016-04-07 12:01:39Z maronga
77! Renamed radiation calls.
78! Renamed canopy model calls.
79!
80! 1822 2016-04-07 07:49:42Z hoffmann
81! icloud_scheme replaced by microphysics_*
82!
83! 1817 2016-04-06 15:44:20Z maronga
84! Renamed lsm calls.
85!
86! 1815 2016-04-06 13:49:59Z raasch
87! zero-settings for velocities inside topography re-activated (was deactivated
88! in r1762)
89!
90! 1788 2016-03-10 11:01:04Z maronga
91! Added z0q.
92! Syntax layout improved.
93!
94! 1783 2016-03-06 18:36:17Z raasch
95! netcdf module name changed + related changes
96!
97! 1764 2016-02-28 12:45:19Z raasch
98! bugfix: increase size of volume_flow_area_l and volume_flow_initial_l by 1
99!
100! 1762 2016-02-25 12:31:13Z hellstea
101! Introduction of nested domain feature
102!
103! 1738 2015-12-18 13:56:05Z raasch
104! calculate mean surface level height for each statistic region
105!
106! 1734 2015-12-02 12:17:12Z raasch
107! no initial disturbances in case that the disturbance energy limit has been
108! set zero
109!
110! 1707 2015-11-02 15:24:52Z maronga
111! Bugfix: transfer of Richardson number from 1D model to Obukhov length caused
112! devision by zero in neutral stratification
113!
114! 1691 2015-10-26 16:17:44Z maronga
115! Call to init_surface_layer added. rif is replaced by ol and zeta.
116!
117! 1682 2015-10-07 23:56:08Z knoop
118! Code annotations made doxygen readable
119!
120! 1615 2015-07-08 18:49:19Z suehring
121! Enable turbulent inflow for passive_scalar and humidity
122!
123! 1585 2015-04-30 07:05:52Z maronga
124! Initialization of radiation code is now done after LSM initializtion
125!
126! 1575 2015-03-27 09:56:27Z raasch
127! adjustments for psolver-queries
128!
129! 1551 2015-03-03 14:18:16Z maronga
130! Allocation of land surface arrays is now done in the subroutine lsm_init_arrays,
131! which is part of land_surface_model.
132!
133! 1507 2014-12-10 12:14:18Z suehring
134! Bugfix: set horizontal velocity components to zero inside topography
135!
136! 1496 2014-12-02 17:25:50Z maronga
137! Added initialization of the land surface and radiation schemes
138!
139! 1484 2014-10-21 10:53:05Z kanani
140! Changes due to new module structure of the plant canopy model:
141! canopy-related initialization (e.g. lad and canopy_heat_flux) moved to new
142! subroutine init_plant_canopy within the module plant_canopy_model_mod,
143! call of subroutine init_plant_canopy added.
144!
145! 1431 2014-07-15 14:47:17Z suehring
146! var_d added, in order to normalize spectra.
147!
148! 1429 2014-07-15 12:53:45Z knoop
149! Ensemble run capability added to parallel random number generator
150!
151! 1411 2014-05-16 18:01:51Z suehring
152! Initial horizontal velocity profiles were not set to zero at the first vertical
153! grid level in case of non-cyclic lateral boundary conditions.
154!
155! 1406 2014-05-16 13:47:01Z raasch
156! bugfix: setting of initial velocities at k=1 to zero not in case of a
157! no-slip boundary condition for uv
158!
159! 1402 2014-05-09 14:25:13Z raasch
160! location messages modified
161!
162! 1400 2014-05-09 14:03:54Z knoop
163! Parallel random number generator added
164!
165! 1384 2014-05-02 14:31:06Z raasch
166! location messages added
167!
168! 1361 2014-04-16 15:17:48Z hoffmann
169! tend_* removed
170! Bugfix: w_subs is not allocated anymore if it is already allocated
171!
172! 1359 2014-04-11 17:15:14Z hoffmann
173! module lpm_init_mod added to use statements, because lpm_init has become a
174! module
175!
176! 1353 2014-04-08 15:21:23Z heinze
177! REAL constants provided with KIND-attribute
178!
179! 1340 2014-03-25 19:45:13Z kanani
180! REAL constants defined as wp-kind
181!
182! 1322 2014-03-20 16:38:49Z raasch
183! REAL constants defined as wp-kind
184! module interfaces removed
185!
186! 1320 2014-03-20 08:40:49Z raasch
187! ONLY-attribute added to USE-statements,
188! kind-parameters added to all INTEGER and REAL declaration statements,
189! kinds are defined in new module kinds,
190! revision history before 2012 removed,
191! comment fields (!:) to be used for variable explanations added to
192! all variable declaration statements
193!
194! 1316 2014-03-17 07:44:59Z heinze
195! Bugfix: allocation of w_subs
196!
197! 1299 2014-03-06 13:15:21Z heinze
198! Allocate w_subs due to extension of large scale subsidence in combination
199! with large scale forcing data (LSF_DATA)
200!
201! 1241 2013-10-30 11:36:58Z heinze
202! Overwrite initial profiles in case of nudging
203! Inititialize shf and qsws in case of large_scale_forcing
204!
205! 1221 2013-09-10 08:59:13Z raasch
206! +rflags_s_inner in copyin statement, use copyin for most arrays instead of
207! copy
208!
209! 1212 2013-08-15 08:46:27Z raasch
210! array tri is allocated and included in data copy statement
211!
212! 1195 2013-07-01 12:27:57Z heinze
213! Bugfix: move allocation of ref_state to parin.f90 and read_var_list.f90
214!
215! 1179 2013-06-14 05:57:58Z raasch
216! allocate and set ref_state to be used in buoyancy terms
217!
218! 1171 2013-05-30 11:27:45Z raasch
219! diss array is allocated with full size if accelerator boards are used
220!
221! 1159 2013-05-21 11:58:22Z fricke
222! -bc_lr_dirneu, bc_lr_neudir, bc_ns_dirneu, bc_ns_neudir
223!
224! 1153 2013-05-10 14:33:08Z raasch
225! diss array is allocated with dummy elements even if it is not needed
226! (required by PGI 13.4 / CUDA 5.0)
227!
228! 1115 2013-03-26 18:16:16Z hoffmann
229! unused variables removed
230!
231! 1113 2013-03-10 02:48:14Z raasch
232! openACC directive modified
233!
234! 1111 2013-03-08 23:54:10Z raasch
235! openACC directives added for pres
236! array diss allocated only if required
237!
238! 1092 2013-02-02 11:24:22Z raasch
239! unused variables removed
240!
241! 1065 2012-11-22 17:42:36Z hoffmann
242! allocation of diss (dissipation rate) in case of turbulence = .TRUE. added
243!
244! 1053 2012-11-13 17:11:03Z hoffmann
245! allocation and initialisation of necessary data arrays for the two-moment
246! cloud physics scheme the two new prognostic equations (nr, qr):
247! +dr, lambda_r, mu_r, sed_*, xr, *s, *sws, *swst, *, *_p, t*_m, *_1, *_2, *_3,
248! +tend_*, prr
249!
250! 1036 2012-10-22 13:43:42Z raasch
251! code put under GPL (PALM 3.9)
252!
253! 1032 2012-10-21 13:03:21Z letzel
254! save memory by not allocating pt_2 in case of neutral = .T.
255!
256! 1025 2012-10-07 16:04:41Z letzel
257! bugfix: swap indices of mask for ghost boundaries
258!
259! 1015 2012-09-27 09:23:24Z raasch
260! mask is set to zero for ghost boundaries
261!
262! 1010 2012-09-20 07:59:54Z raasch
263! cpp switch __nopointer added for pointer free version
264!
265! 1003 2012-09-14 14:35:53Z raasch
266! nxra,nyna, nzta replaced ny nxr, nyn, nzt
267!
268! 1001 2012-09-13 14:08:46Z raasch
269! all actions concerning leapfrog scheme removed
270!
271! 996 2012-09-07 10:41:47Z raasch
272! little reformatting
273!
274! 978 2012-08-09 08:28:32Z fricke
275! outflow damping layer removed
276! roughness length for scalar quantites z0h added
277! damping zone for the potential temperatur in case of non-cyclic lateral
278! boundaries added
279! initialization of ptdf_x, ptdf_y
280! initialization of c_u_m, c_u_m_l, c_v_m, c_v_m_l, c_w_m, c_w_m_l
281!
282! 849 2012-03-15 10:35:09Z raasch
283! init_particles renamed lpm_init
284!
285! 825 2012-02-19 03:03:44Z raasch
286! wang_collision_kernel renamed wang_kernel
287!
288! Revision 1.1  1998/03/09 16:22:22  raasch
289! Initial revision
290!
291!
292! Description:
293! ------------
294!> Allocation of arrays and initialization of the 3D model via
295!> a) pre-run the 1D model
296!> or
297!> b) pre-set constant linear profiles
298!> or
299!> c) read values of a previous run
300!------------------------------------------------------------------------------!
301 SUBROUTINE init_3d_model
302 
303
304    USE advec_ws
305
306    USE arrays_3d
307
308    USE constants,                                                             &
309        ONLY:  pi
310   
311    USE control_parameters
312   
313    USE flight_mod,                                                            &
314        ONLY:  flight_init
315   
316    USE grid_variables,                                                        &
317        ONLY:  dx, dy
318   
319    USE indices
320
321    USE lpm_init_mod,                                                          &
322        ONLY:  lpm_init
323   
324    USE kinds
325
326    USE land_surface_model_mod,                                                &
327        ONLY:  lsm_init, lsm_init_arrays, land_surface
328 
329    USE ls_forcing_mod
330
331    USE microphysics_mod,                                                      &
332        ONLY:  collision_turbulence, microphysics_init
333
334    USE model_1d,                                                              &
335        ONLY:  e1d, kh1d, km1d, l1d, rif1d, u1d, us1d, usws1d, v1d, vsws1d
336   
337    USE netcdf_interface,                                                      &
338        ONLY:  dots_max, dots_num
339   
340    USE particle_attributes,                                                   &
341        ONLY:  particle_advection, use_sgs_for_particles, wang_kernel
342   
343    USE pegrid
344   
345    USE plant_canopy_model_mod,                                                &
346        ONLY:  pcm_init, plant_canopy
347
348    USE radiation_model_mod,                                                   &
349        ONLY:  radiation_init, radiation
350   
351    USE random_function_mod
352   
353    USE random_generator_parallel,                                             &
354        ONLY:  random_number_parallel, random_seed_parallel, random_dummy,     &
355               id_random_array, seq_random_array
356   
357    USE statistics,                                                            &
358        ONLY:  hom, hom_sum, mean_surface_level_height, pr_palm, rmask,        &
359               statistic_regions, sums, sums_divnew_l, sums_divold_l, sums_l,  &
360               sums_l_l, sums_up_fraction_l, sums_wsts_bc_l, ts_value,         &
361               weight_pres, weight_substep
362 
363    USE surface_layer_fluxes_mod,                                              &
364        ONLY:  init_surface_layer_fluxes
365   
366    USE transpose_indices
367
368    USE urban_surface_mod,                                                     &
369        ONLY:  urban_surface, usm_init_urban_surface
370
371    USE wind_turbine_model_mod,                                                &
372        ONLY:  wtm_init, wtm_init_arrays, wind_turbine
373
374    IMPLICIT NONE
375
376    INTEGER(iwp) ::  i             !<
377    INTEGER(iwp) ::  ind_array(1)  !<
378    INTEGER(iwp) ::  j             !<
379    INTEGER(iwp) ::  k             !<
380    INTEGER(iwp) ::  sr            !<
381
382    INTEGER(iwp), DIMENSION(:), ALLOCATABLE   ::  ngp_2dh_l  !<
383
384    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  ngp_2dh_outer_l    !<
385    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  ngp_2dh_s_inner_l  !<
386
387    REAL(wp), DIMENSION(1:3) ::  volume_flow_area_l     !<
388    REAL(wp), DIMENSION(1:3) ::  volume_flow_initial_l  !<
389
390    REAL(wp), DIMENSION(:), ALLOCATABLE ::  mean_surface_level_height_l    !<
391    REAL(wp), DIMENSION(:), ALLOCATABLE ::  ngp_3d_inner_l    !<
392    REAL(wp), DIMENSION(:), ALLOCATABLE ::  ngp_3d_inner_tmp  !<
393
394
395    CALL location_message( 'allocating arrays', .FALSE. )
396!
397!-- Allocate arrays
398    ALLOCATE( mean_surface_level_height(0:statistic_regions),                  &
399              mean_surface_level_height_l(0:statistic_regions),                &
400              ngp_2dh(0:statistic_regions), ngp_2dh_l(0:statistic_regions),    &
401              ngp_3d(0:statistic_regions),                                     &
402              ngp_3d_inner(0:statistic_regions),                               &
403              ngp_3d_inner_l(0:statistic_regions),                             &
404              ngp_3d_inner_tmp(0:statistic_regions),                           &
405              sums_divnew_l(0:statistic_regions),                              &
406              sums_divold_l(0:statistic_regions) )
407    ALLOCATE( dp_smooth_factor(nzb:nzt), rdf(nzb+1:nzt), rdf_sc(nzb+1:nzt) )
408    ALLOCATE( ngp_2dh_outer(nzb:nzt+1,0:statistic_regions),                    &
409              ngp_2dh_outer_l(nzb:nzt+1,0:statistic_regions),                  &
410              ngp_2dh_s_inner(nzb:nzt+1,0:statistic_regions),                  &
411              ngp_2dh_s_inner_l(nzb:nzt+1,0:statistic_regions),                &
412              rmask(nysg:nyng,nxlg:nxrg,0:statistic_regions),                  &
413              sums(nzb:nzt+1,pr_palm+max_pr_user),                             &
414              sums_l(nzb:nzt+1,pr_palm+max_pr_user,0:threads_per_task-1),      &
415              sums_l_l(nzb:nzt+1,0:statistic_regions,0:threads_per_task-1),    &
416              sums_up_fraction_l(10,3,0:statistic_regions),                    &
417              sums_wsts_bc_l(nzb:nzt+1,0:statistic_regions),                   &
418              ts_value(dots_max,0:statistic_regions) )
419    ALLOCATE( ptdf_x(nxlg:nxrg), ptdf_y(nysg:nyng) )
420
421    ALLOCATE( ol(nysg:nyng,nxlg:nxrg), shf(nysg:nyng,nxlg:nxrg),               &
422              ts(nysg:nyng,nxlg:nxrg), tswst(nysg:nyng,nxlg:nxrg),             &
423              us(nysg:nyng,nxlg:nxrg), usws(nysg:nyng,nxlg:nxrg),              &
424              uswst(nysg:nyng,nxlg:nxrg), vsws(nysg:nyng,nxlg:nxrg),           &
425              vswst(nysg:nyng,nxlg:nxrg), z0(nysg:nyng,nxlg:nxrg),             &
426              z0h(nysg:nyng,nxlg:nxrg), z0q(nysg:nyng,nxlg:nxrg) )
427
428    ALLOCATE( d(nzb+1:nzt,nys:nyn,nxl:nxr),                                    &
429              kh(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                               &
430              km(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                               &
431              p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                &
432              tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
433
434#if defined( __nopointer )
435    ALLOCATE( e(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                &
436              e_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
437              pt(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                               &
438              pt_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
439              u(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                &
440              u_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
441              v(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                &
442              v_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
443              w(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                &
444              w_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
445              te_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
446              tpt_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                            &
447              tu_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
448              tv_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
449              tw_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
450#else
451    ALLOCATE( e_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
452              e_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
453              e_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
454              pt_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
455              pt_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
456              u_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
457              u_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
458              u_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
459              v_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
460              v_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
461              v_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
462              w_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
463              w_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
464              w_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
465    IF (  .NOT.  neutral )  THEN
466       ALLOCATE( pt_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
467    ENDIF
468#endif
469
470!
471!-- Following array is required for perturbation pressure within the iterative
472!-- pressure solvers. For the multistep schemes (Runge-Kutta), array p holds
473!-- the weighted average of the substeps and cannot be used in the Poisson
474!-- solver.
475    IF ( psolver == 'sor' )  THEN
476       ALLOCATE( p_loc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
477    ELSEIF ( psolver(1:9) == 'multigrid' )  THEN
478!
479!--    For performance reasons, multigrid is using one ghost layer only
480       ALLOCATE( p_loc(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
481    ENDIF
482
483!
484!-- Array for storing constant coeffficients of the tridiagonal solver
485    IF ( psolver == 'poisfft' )  THEN
486       ALLOCATE( tri(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1,2) )
487       ALLOCATE( tric(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1) )
488    ENDIF
489
490    IF ( humidity )  THEN
491!
492!--    2D-humidity
493       ALLOCATE ( qs(nysg:nyng,nxlg:nxrg),                                     &
494                  qsws(nysg:nyng,nxlg:nxrg),                                   &
495                  qswst(nysg:nyng,nxlg:nxrg) )
496
497!
498!--    3D-humidity
499#if defined( __nopointer )
500       ALLOCATE( q(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
501                 q_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
502                 tq_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
503#else
504       ALLOCATE( q_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
505                 q_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
506                 q_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
507#endif
508
509!
510!--    3D-arrays needed for humidity
511       IF ( humidity )  THEN
512#if defined( __nopointer )
513          ALLOCATE( vpt(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
514#else
515          ALLOCATE( vpt_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
516#endif
517
518          IF ( cloud_physics )  THEN
519
520!
521!--          Liquid water content
522#if defined( __nopointer )
523             ALLOCATE ( ql(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
524#else
525             ALLOCATE ( ql_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
526#endif
527!
528!--          Precipitation amount and rate (only needed if output is switched)
529             ALLOCATE( precipitation_amount(nysg:nyng,nxlg:nxrg),              &
530                       precipitation_rate(nysg:nyng,nxlg:nxrg) )
531
532!
533!--          3D-cloud water content
534#if defined( __nopointer )
535             ALLOCATE( qc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
536#else
537             ALLOCATE( qc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
538#endif
539!
540!--          3d-precipitation rate
541             ALLOCATE( prr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
542
543             IF ( microphysics_seifert )  THEN
544!
545!--             2D-rain water content and rain drop concentration arrays
546                ALLOCATE ( qrs(nysg:nyng,nxlg:nxrg),                        &
547                           qrsws(nysg:nyng,nxlg:nxrg),                      &
548                           qrswst(nysg:nyng,nxlg:nxrg),                     &
549                           nrs(nysg:nyng,nxlg:nxrg),                        &
550                           nrsws(nysg:nyng,nxlg:nxrg),                      &
551                           nrswst(nysg:nyng,nxlg:nxrg) )
552!
553!--             3D-rain water content, rain drop concentration arrays
554#if defined( __nopointer )
555                ALLOCATE( nr(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                &
556                          nr_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
557                          qr(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                &
558                          qr_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
559                          tnr_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg),             &
560                          tqr_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
561#else
562                ALLOCATE( nr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
563                          nr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
564                          nr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
565                          qr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
566                          qr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
567                          qr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
568#endif
569             ENDIF
570
571          ENDIF
572
573          IF ( cloud_droplets )  THEN
574!
575!--          Liquid water content, change in liquid water content
576#if defined( __nopointer )
577             ALLOCATE ( ql(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                     &
578                        ql_c(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
579#else
580             ALLOCATE ( ql_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                   &
581                        ql_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
582#endif
583!
584!--          Real volume of particles (with weighting), volume of particles
585             ALLOCATE ( ql_v(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                   &
586                        ql_vp(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
587          ENDIF
588
589       ENDIF
590
591    ENDIF
592   
593   
594    IF ( passive_scalar )  THEN
595!
596!--    2D-scalar arrays
597       ALLOCATE ( ss(nysg:nyng,nxlg:nxrg),                                     &
598                  ssws(nysg:nyng,nxlg:nxrg),                                   &
599                  sswst(nysg:nyng,nxlg:nxrg) )
600
601!
602!--    3D scalar arrays
603#if defined( __nopointer )
604       ALLOCATE( s(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
605                 s_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
606                 ts_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
607#else
608       ALLOCATE( s_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
609                 s_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
610                 s_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
611#endif
612    ENDIF
613
614    IF ( ocean )  THEN
615       ALLOCATE( saswsb(nysg:nyng,nxlg:nxrg),                                  &
616                 saswst(nysg:nyng,nxlg:nxrg) )
617#if defined( __nopointer )
618       ALLOCATE( prho(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                          &
619                 rho(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
620                 sa(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                            &
621                 sa_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                          &
622                 tsa_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
623#else
624       ALLOCATE( prho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                        &
625                 rho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                         &
626                 sa_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                          &
627                 sa_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                          &
628                 sa_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
629       prho => prho_1
630       rho  => rho_1  ! routines calc_mean_profile and diffusion_e require
631                      ! density to be apointer
632#endif
633       IF ( humidity_remote )  THEN
634          ALLOCATE( qswst_remote(nysg:nyng,nxlg:nxrg))
635          qswst_remote = 0.0_wp
636       ENDIF
637    ENDIF
638
639!
640!-- 3D-array for storing the dissipation, needed for calculating the sgs
641!-- particle velocities
642    IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.  collision_turbulence  &
643         .OR.  num_acc_per_node > 0 )  THEN
644       ALLOCATE( diss(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
645    ENDIF
646
647!
648!-- 1D-array for large scale subsidence velocity
649    IF ( .NOT. ALLOCATED( w_subs ) )  THEN
650       ALLOCATE ( w_subs(nzb:nzt+1) )
651       w_subs = 0.0_wp
652    ENDIF
653
654!
655!-- ID-array and state-space-array for the parallel random number generator
656    IF ( random_generator == 'random-parallel' )  THEN
657       ALLOCATE ( seq_random_array(5,nysg:nyng,nxlg:nxrg) )
658       ALLOCATE ( id_random_array(0:ny,0:nx) )
659       seq_random_array = 0
660       id_random_array  = 0
661    ENDIF
662   
663!
664!-- 4D-array for storing the Rif-values at vertical walls
665    IF ( topography /= 'flat' )  THEN
666       ALLOCATE( rif_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg,1:4) )
667       rif_wall = 0.0_wp
668    ENDIF
669
670!
671!-- Arrays to store velocity data from t-dt and the phase speeds which
672!-- are needed for radiation boundary conditions
673    IF ( outflow_l )  THEN
674       ALLOCATE( u_m_l(nzb:nzt+1,nysg:nyng,1:2),                               &
675                 v_m_l(nzb:nzt+1,nysg:nyng,0:1),                               &
676                 w_m_l(nzb:nzt+1,nysg:nyng,0:1) )
677    ENDIF
678    IF ( outflow_r )  THEN
679       ALLOCATE( u_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx),                           &
680                 v_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx),                           &
681                 w_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx) )
682    ENDIF
683    IF ( outflow_l  .OR.  outflow_r )  THEN
684       ALLOCATE( c_u(nzb:nzt+1,nysg:nyng), c_v(nzb:nzt+1,nysg:nyng),           &
685                 c_w(nzb:nzt+1,nysg:nyng) )
686    ENDIF
687    IF ( outflow_s )  THEN
688       ALLOCATE( u_m_s(nzb:nzt+1,0:1,nxlg:nxrg),                               &
689                 v_m_s(nzb:nzt+1,1:2,nxlg:nxrg),                               &
690                 w_m_s(nzb:nzt+1,0:1,nxlg:nxrg) )
691    ENDIF
692    IF ( outflow_n )  THEN
693       ALLOCATE( u_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg),                           &
694                 v_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg),                           &
695                 w_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg) )
696    ENDIF
697    IF ( outflow_s  .OR.  outflow_n )  THEN
698       ALLOCATE( c_u(nzb:nzt+1,nxlg:nxrg), c_v(nzb:nzt+1,nxlg:nxrg),           &
699                 c_w(nzb:nzt+1,nxlg:nxrg) )
700    ENDIF
701    IF ( outflow_l  .OR.  outflow_r  .OR.  outflow_s  .OR.  outflow_n )  THEN
702       ALLOCATE( c_u_m_l(nzb:nzt+1), c_v_m_l(nzb:nzt+1), c_w_m_l(nzb:nzt+1) )                   
703       ALLOCATE( c_u_m(nzb:nzt+1), c_v_m(nzb:nzt+1), c_w_m(nzb:nzt+1) )
704    ENDIF
705
706
707#if ! defined( __nopointer )
708!
709!-- Initial assignment of the pointers
710    e  => e_1;   e_p  => e_2;   te_m  => e_3
711    IF ( .NOT. neutral )  THEN
712       pt => pt_1;  pt_p => pt_2;  tpt_m => pt_3
713    ELSE
714       pt => pt_1;  pt_p => pt_1;  tpt_m => pt_3
715    ENDIF
716    u  => u_1;   u_p  => u_2;   tu_m  => u_3
717    v  => v_1;   v_p  => v_2;   tv_m  => v_3
718    w  => w_1;   w_p  => w_2;   tw_m  => w_3
719
720    IF ( humidity )  THEN
721       q => q_1;  q_p => q_2;  tq_m => q_3
722       IF ( humidity )  THEN
723          vpt  => vpt_1   
724          IF ( cloud_physics )  THEN
725             ql => ql_1
726             qc => qc_1
727             IF ( microphysics_seifert )  THEN
728                qr => qr_1;  qr_p  => qr_2;  tqr_m  => qr_3
729                nr => nr_1;  nr_p  => nr_2;  tnr_m  => nr_3
730             ENDIF
731          ENDIF
732       ENDIF
733       IF ( cloud_droplets )  THEN
734          ql   => ql_1
735          ql_c => ql_2
736       ENDIF
737    ENDIF
738   
739    IF ( passive_scalar )  THEN
740       s => s_1;  s_p => s_2;  ts_m => s_3
741    ENDIF   
742
743    IF ( ocean )  THEN
744       sa => sa_1;  sa_p => sa_2;  tsa_m => sa_3
745    ENDIF
746#endif
747
748!
749!-- Allocate land surface model arrays
750    IF ( land_surface )  THEN
751       CALL lsm_init_arrays
752    ENDIF
753
754!
755!-- Allocate wind turbine model arrays
756    IF ( wind_turbine )  THEN
757       CALL wtm_init_arrays
758    ENDIF
759   
760!
761!-- Initialize virtual flight measurements
762    IF ( virtual_flight )  THEN
763       CALL flight_init
764    ENDIF
765
766!
767!-- Allocate arrays containing the RK coefficient for calculation of
768!-- perturbation pressure and turbulent fluxes. At this point values are
769!-- set for pressure calculation during initialization (where no timestep
770!-- is done). Further below the values needed within the timestep scheme
771!-- will be set.
772    ALLOCATE( weight_substep(1:intermediate_timestep_count_max),               &
773              weight_pres(1:intermediate_timestep_count_max) )
774    weight_substep = 1.0_wp
775    weight_pres    = 1.0_wp
776    intermediate_timestep_count = 0  ! needed when simulated_time = 0.0
777       
778    CALL location_message( 'finished', .TRUE. )
779
780!
781!-- Initialize local summation arrays for routine flow_statistics.
782!-- This is necessary because they may not yet have been initialized when they
783!-- are called from flow_statistics (or - depending on the chosen model run -
784!-- are never initialized)
785    sums_divnew_l      = 0.0_wp
786    sums_divold_l      = 0.0_wp
787    sums_l_l           = 0.0_wp
788    sums_up_fraction_l = 0.0_wp
789    sums_wsts_bc_l     = 0.0_wp
790
791
792!
793!-- Initialize model variables
794    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
795         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
796!
797!--    First model run of a possible job queue.
798!--    Initial profiles of the variables must be computes.
799       IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
800
801          CALL location_message( 'initializing with 1D model profiles', .FALSE. )
802!
803!--       Use solutions of the 1D model as initial profiles,
804!--       start 1D model
805          CALL init_1d_model
806!
807!--       Transfer initial profiles to the arrays of the 3D model
808          DO  i = nxlg, nxrg
809             DO  j = nysg, nyng
810                e(:,j,i)  = e1d
811                kh(:,j,i) = kh1d
812                km(:,j,i) = km1d
813                pt(:,j,i) = pt_init
814                u(:,j,i)  = u1d
815                v(:,j,i)  = v1d
816             ENDDO
817          ENDDO
818
819          IF ( humidity )  THEN
820             DO  i = nxlg, nxrg
821                DO  j = nysg, nyng
822                   q(:,j,i) = q_init
823                ENDDO
824             ENDDO
825             IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
826                DO  i = nxlg, nxrg
827                   DO  j = nysg, nyng
828                      qr(:,j,i) = 0.0_wp
829                      nr(:,j,i) = 0.0_wp
830                   ENDDO
831                ENDDO
832
833             ENDIF
834          ENDIF
835          IF ( passive_scalar )  THEN
836             DO  i = nxlg, nxrg
837                DO  j = nysg, nyng
838                   s(:,j,i) = s_init
839                ENDDO
840             ENDDO   
841          ENDIF
842
843          IF ( .NOT. constant_diffusion )  THEN
844             DO  i = nxlg, nxrg
845                DO  j = nysg, nyng
846                   e(:,j,i)  = e1d
847                ENDDO
848             ENDDO
849!
850!--          Store initial profiles for output purposes etc.
851             hom(:,1,25,:) = SPREAD( l1d, 2, statistic_regions+1 )
852
853             IF ( constant_flux_layer )  THEN
854                ol   = ( zu(nzb+1) - zw(nzb) ) / ( rif1d(nzb+1) + 1.0E-20_wp )
855                ts   = 0.0_wp  ! could actually be computed more accurately in the
856                               ! 1D model. Update when opportunity arises.
857                us   = us1d
858                usws = usws1d
859                vsws = vsws1d
860             ELSE
861                ts   = 0.0_wp  ! must be set, because used in
862                ol   = ( zu(nzb+1) - zw(nzb) ) / zeta_min  ! flowste
863                us   = 0.0_wp
864                usws = 0.0_wp
865                vsws = 0.0_wp
866             ENDIF
867
868          ELSE
869             e    = 0.0_wp  ! must be set, because used in
870             ol   = ( zu(nzb+1) - zw(nzb) ) / zeta_min  ! flowste
871             ts   = 0.0_wp
872             us   = 0.0_wp
873             usws = 0.0_wp
874             vsws = 0.0_wp
875          ENDIF
876          uswst = top_momentumflux_u
877          vswst = top_momentumflux_v
878
879!
880!--       In every case qs = 0.0 (see also pt)
881!--       This could actually be computed more accurately in the 1D model.
882!--       Update when opportunity arises!
883          IF ( humidity )  THEN
884             qs = 0.0_wp
885             IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
886                qrs = 0.0_wp
887                nrs = 0.0_wp
888             ENDIF
889          ENDIF
890!
891!--       Initialize scaling parameter for passive scalar
892          IF ( passive_scalar ) ss = 0.0_wp         
893
894!
895!--       Inside buildings set velocities back to zero
896          IF ( topography /= 'flat' )  THEN
897             DO  i = nxlg, nxrg
898                DO  j = nysg, nyng
899                   u(nzb:nzb_u_inner(j,i),j,i) = 0.0_wp
900                   v(nzb:nzb_v_inner(j,i),j,i) = 0.0_wp
901                ENDDO
902             ENDDO
903             
904!
905!--          WARNING: The extra boundary conditions set after running the
906!--          -------  1D model impose an error on the divergence one layer
907!--                   below the topography; need to correct later
908!--          ATTENTION: Provisional correction for Piacsek & Williams
909!--          ---------  advection scheme: keep u and v zero one layer below
910!--                     the topography.
911             IF ( ibc_uv_b == 1 )  THEN
912!
913!--             Neumann condition
914                DO  i = nxl-1, nxr+1
915                   DO  j = nys-1, nyn+1
916                      IF ( nzb_u_inner(j,i) == 0 ) u(0,j,i) = u(1,j,i)
917                      IF ( nzb_v_inner(j,i) == 0 ) v(0,j,i) = v(1,j,i)
918                   ENDDO
919                ENDDO
920
921             ENDIF
922
923          ENDIF
924
925          CALL location_message( 'finished', .TRUE. )
926
927       ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 )    &
928       THEN
929
930          CALL location_message( 'initializing with constant profiles', .FALSE. )
931!
932!--       Overwrite initial profiles in case of nudging
933          IF ( nudging )  THEN
934             pt_init = ptnudge(:,1)
935             u_init  = unudge(:,1)
936             v_init  = vnudge(:,1)
937             IF ( humidity  )  THEN ! is passive_scalar correct???
938                q_init = qnudge(:,1)
939             ENDIF
940
941             WRITE( message_string, * ) 'Initial profiles of u, v and ',       &
942                 'scalars from NUDGING_DATA are used.'
943             CALL message( 'init_3d_model', 'PA0370', 0, 0, 0, 6, 0 )
944          ENDIF
945
946!
947!--       Use constructed initial profiles (velocity constant with height,
948!--       temperature profile with constant gradient)
949          DO  i = nxlg, nxrg
950             DO  j = nysg, nyng
951                pt(:,j,i) = pt_init
952                u(:,j,i)  = u_init
953                v(:,j,i)  = v_init
954             ENDDO
955          ENDDO
956
957!
958!--       Set initial horizontal velocities at the lowest computational grid
959!--       levels to zero in order to avoid too small time steps caused by the
960!--       diffusion limit in the initial phase of a run (at k=1, dz/2 occurs
961!--       in the limiting formula!).
962          IF ( ibc_uv_b /= 1 )  THEN
963             DO  i = nxlg, nxrg
964                DO  j = nysg, nyng
965                   u(nzb:nzb_u_inner(j,i)+1,j,i) = 0.0_wp
966                   v(nzb:nzb_v_inner(j,i)+1,j,i) = 0.0_wp
967                ENDDO
968             ENDDO
969          ENDIF
970
971          IF ( humidity )  THEN
972             DO  i = nxlg, nxrg
973                DO  j = nysg, nyng
974                   q(:,j,i) = q_init
975                ENDDO
976             ENDDO
977             IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
978
979                DO  i = nxlg, nxrg
980                   DO  j = nysg, nyng
981                      qr(:,j,i) = 0.0_wp
982                      nr(:,j,i) = 0.0_wp
983                   ENDDO
984                ENDDO
985
986             ENDIF
987          ENDIF
988         
989          IF ( passive_scalar )  THEN
990             DO  i = nxlg, nxrg
991                DO  j = nysg, nyng
992                   s(:,j,i) = s_init
993                ENDDO
994             ENDDO
995          ENDIF
996
997          IF ( ocean )  THEN
998             DO  i = nxlg, nxrg
999                DO  j = nysg, nyng
1000                   sa(:,j,i) = sa_init
1001                ENDDO
1002             ENDDO
1003          ENDIF
1004         
1005          IF ( constant_diffusion )  THEN
1006             km   = km_constant
1007             kh   = km / prandtl_number
1008             e    = 0.0_wp
1009          ELSEIF ( e_init > 0.0_wp )  THEN
1010             DO  k = nzb+1, nzt
1011                km(k,:,:) = 0.1_wp * l_grid(k) * SQRT( e_init )
1012             ENDDO
1013             km(nzb,:,:)   = km(nzb+1,:,:)
1014             km(nzt+1,:,:) = km(nzt,:,:)
1015             kh   = km / prandtl_number
1016             e    = e_init
1017          ELSE
1018             IF ( .NOT. ocean )  THEN
1019                kh   = 0.01_wp   ! there must exist an initial diffusion, because
1020                km   = 0.01_wp   ! otherwise no TKE would be produced by the
1021                              ! production terms, as long as not yet
1022                              ! e = (u*/cm)**2 at k=nzb+1
1023             ELSE
1024                kh   = 0.00001_wp
1025                km   = 0.00001_wp
1026             ENDIF
1027             e    = 0.0_wp
1028          ENDIF
1029          ol    = ( zu(nzb+1) - zw(nzb) ) / zeta_min
1030          ts    = 0.0_wp
1031!
1032!--       Very small number is required for calculation of Obukhov length
1033!--       at first timestep     
1034          us    = 1E-30_wp 
1035          usws  = 0.0_wp
1036          uswst = top_momentumflux_u
1037          vsws  = 0.0_wp
1038          vswst = top_momentumflux_v
1039          IF ( humidity )       qs = 0.0_wp
1040          IF ( passive_scalar ) ss = 0.0_wp
1041
1042!
1043!--       Compute initial temperature field and other constants used in case
1044!--       of a sloping surface
1045          IF ( sloping_surface )  CALL init_slope
1046
1047          CALL location_message( 'finished', .TRUE. )
1048
1049       ELSEIF ( INDEX(initializing_actions, 'by_user') /= 0 )                  &
1050       THEN
1051
1052          CALL location_message( 'initializing by user', .FALSE. )
1053!
1054!--       Initialization will completely be done by the user
1055          CALL user_init_3d_model
1056
1057          CALL location_message( 'finished', .TRUE. )
1058
1059       ENDIF
1060
1061       CALL location_message( 'initializing statistics, boundary conditions, etc.', &
1062                              .FALSE. )
1063
1064!
1065!--    Bottom boundary
1066       IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2  )  THEN
1067          u(nzb,:,:) = 0.0_wp
1068          v(nzb,:,:) = 0.0_wp
1069       ENDIF
1070
1071!
1072!--    Apply channel flow boundary condition
1073       IF ( TRIM( bc_uv_t ) == 'dirichlet_0' )  THEN
1074          u(nzt+1,:,:) = 0.0_wp
1075          v(nzt+1,:,:) = 0.0_wp
1076       ENDIF
1077
1078!
1079!--    Calculate virtual potential temperature
1080       IF ( humidity )  vpt = pt * ( 1.0_wp + 0.61_wp * q )
1081
1082!
1083!--    Store initial profiles for output purposes etc.
1084       hom(:,1,5,:) = SPREAD( u(:,nys,nxl), 2, statistic_regions+1 )
1085       hom(:,1,6,:) = SPREAD( v(:,nys,nxl), 2, statistic_regions+1 )
1086       IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2)  THEN
1087          hom(nzb,1,5,:) = 0.0_wp
1088          hom(nzb,1,6,:) = 0.0_wp
1089       ENDIF
1090       hom(:,1,7,:)  = SPREAD( pt(:,nys,nxl), 2, statistic_regions+1 )
1091       hom(:,1,23,:) = SPREAD( km(:,nys,nxl), 2, statistic_regions+1 )
1092       hom(:,1,24,:) = SPREAD( kh(:,nys,nxl), 2, statistic_regions+1 )
1093
1094       IF ( ocean )  THEN
1095!
1096!--       Store initial salinity profile
1097          hom(:,1,26,:)  = SPREAD( sa(:,nys,nxl), 2, statistic_regions+1 )
1098       ENDIF
1099
1100       IF ( humidity )  THEN
1101!
1102!--       Store initial profile of total water content, virtual potential
1103!--       temperature
1104          hom(:,1,26,:) = SPREAD(   q(:,nys,nxl), 2, statistic_regions+1 )
1105          hom(:,1,29,:) = SPREAD( vpt(:,nys,nxl), 2, statistic_regions+1 )
1106          IF ( cloud_physics  .OR.  cloud_droplets ) THEN
1107!
1108!--          Store initial profile of specific humidity and potential
1109!--          temperature
1110             hom(:,1,27,:) = SPREAD(  q(:,nys,nxl), 2, statistic_regions+1 )
1111             hom(:,1,28,:) = SPREAD( pt(:,nys,nxl), 2, statistic_regions+1 )
1112          ENDIF
1113       ENDIF
1114
1115       IF ( passive_scalar )  THEN
1116!
1117!--       Store initial scalar profile
1118          hom(:,1,115,:) = SPREAD(  s(:,nys,nxl), 2, statistic_regions+1 )
1119       ENDIF
1120
1121!
1122!--    Initialize the random number generators (from numerical recipes)
1123       CALL random_function_ini
1124       
1125       IF ( random_generator == 'random-parallel' )  THEN
1126!--       Asigning an ID to every vertical gridpoint column
1127!--       dependig on the ensemble run number.
1128          random_dummy=1
1129          DO j=0,ny
1130             DO i=0,nx
1131                id_random_array(j,i) = random_dummy + 1E6                      &
1132                                       * ( ensemble_member_nr - 1000 )
1133                random_dummy = random_dummy + 1
1134             END DO
1135          ENDDO
1136!--       Initializing with random_seed_parallel for every vertical
1137!--       gridpoint column.
1138          random_dummy=0
1139          DO j = nysg, nyng
1140             DO i = nxlg, nxrg
1141                CALL random_seed_parallel (random_sequence=id_random_array(j, i))
1142                CALL random_number_parallel (random_dummy)
1143                CALL random_seed_parallel (get=seq_random_array(:, j, i))
1144             END DO
1145          ENDDO
1146       ENDIF
1147
1148!
1149!--    Initialize fluxes at bottom surface
1150       IF ( use_surface_fluxes )  THEN
1151
1152          IF ( constant_heatflux )  THEN
1153!
1154!--          Heat flux is prescribed
1155             IF ( random_heatflux )  THEN
1156                CALL disturb_heatflux
1157             ELSE
1158                shf = surface_heatflux
1159!
1160!--             Initialize shf with data from external file LSF_DATA
1161                IF ( large_scale_forcing .AND. lsf_surf )  THEN
1162                   CALL ls_forcing_surf ( simulated_time )
1163                ENDIF
1164
1165!
1166!--             Over topography surface_heatflux is replaced by wall_heatflux(0)
1167                IF ( TRIM( topography ) /= 'flat' )  THEN
1168                   DO  i = nxlg, nxrg
1169                      DO  j = nysg, nyng
1170                         IF ( nzb_s_inner(j,i) /= 0 )  THEN
1171                            shf(j,i) = wall_heatflux(0)
1172                         ENDIF
1173                      ENDDO
1174                   ENDDO
1175                ENDIF
1176             ENDIF
1177          ENDIF
1178
1179!
1180!--       Determine the near-surface water flux
1181          IF ( humidity )  THEN
1182             IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
1183                qrsws = 0.0_wp
1184                nrsws = 0.0_wp
1185             ENDIF
1186             IF ( constant_waterflux )  THEN
1187                qsws   = surface_waterflux
1188!
1189!--             Over topography surface_waterflux is replaced by
1190!--             wall_humidityflux(0)
1191                IF ( TRIM( topography ) /= 'flat' )  THEN
1192                   wall_qflux = wall_humidityflux
1193                   DO  i = nxlg, nxrg
1194                      DO  j = nysg, nyng
1195                         IF ( nzb_s_inner(j,i) /= 0 )  THEN
1196                            qsws(j,i) = wall_qflux(0)
1197                         ENDIF
1198                      ENDDO
1199                   ENDDO
1200                ENDIF
1201             ENDIF
1202          ENDIF
1203!
1204!--       Initialize the near-surface scalar flux
1205          IF ( passive_scalar )  THEN
1206             IF ( constant_scalarflux )  THEN
1207                ssws   = surface_scalarflux
1208!
1209!--             Over topography surface_scalarflux is replaced by
1210!--             wall_scalarflux(0)
1211                IF ( TRIM( topography ) /= 'flat' )  THEN
1212                   wall_sflux = wall_scalarflux
1213                   DO  i = nxlg, nxrg
1214                      DO  j = nysg, nyng
1215                         IF ( nzb_s_inner(j,i) /= 0 )  ssws(j,i) = wall_sflux(0)
1216                      ENDDO
1217                   ENDDO
1218                ENDIF
1219             ENDIF
1220          ENDIF   
1221!
1222!--       Initialize near-surface salinity flux
1223          IF ( ocean )  saswsb = bottom_salinityflux
1224
1225       ENDIF
1226
1227!
1228!--    Initialize fluxes at top surface
1229!--    Currently, only the heatflux and salinity flux can be prescribed.
1230!--    The latent flux is zero in this case!
1231       IF ( use_top_fluxes )  THEN
1232!
1233!--       Prescribe to heat flux
1234          IF ( constant_top_heatflux )  tswst = top_heatflux
1235!
1236!--       Prescribe zero latent flux at the top     
1237          IF ( humidity )  THEN
1238             qswst = 0.0_wp
1239             IF ( cloud_physics  .AND.  microphysics_seifert ) THEN
1240                nrswst = 0.0_wp
1241                qrswst = 0.0_wp
1242             ENDIF
1243          ENDIF
1244!
1245!--       Prescribe top scalar flux
1246          IF ( passive_scalar .AND. constant_top_scalarflux )                  &
1247             sswst = top_scalarflux
1248!
1249!--       Prescribe top salinity flux
1250          IF ( ocean .AND. constant_top_salinityflux)                          &
1251             saswst = top_salinityflux
1252!
1253!--       Initialization in case of a coupled model run
1254          IF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
1255             tswst = 0.0_wp
1256          ENDIF
1257
1258       ENDIF
1259
1260!
1261!--    Initialize Prandtl layer quantities
1262       IF ( constant_flux_layer )  THEN
1263
1264          z0 = roughness_length
1265          z0h = z0h_factor * z0
1266          z0q = z0h_factor * z0
1267
1268          IF ( .NOT. constant_heatflux )  THEN 
1269!
1270!--          Surface temperature is prescribed. Here the heat flux cannot be
1271!--          simply estimated, because therefore ol, u* and theta* would have
1272!--          to be computed by iteration. This is why the heat flux is assumed
1273!--          to be zero before the first time step. It approaches its correct
1274!--          value in the course of the first few time steps.
1275             shf   = 0.0_wp
1276          ENDIF
1277
1278          IF ( humidity  )  THEN
1279             IF (  .NOT.  constant_waterflux )  qsws   = 0.0_wp
1280             IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
1281                qrsws = 0.0_wp
1282                nrsws = 0.0_wp
1283             ENDIF
1284          ENDIF
1285          IF ( passive_scalar  .AND.  .NOT.  constant_scalarflux )  ssws = 0.0_wp
1286
1287       ENDIF
1288
1289!
1290!--    Set the reference state to be used in the buoyancy terms (for ocean runs
1291!--    the reference state will be set (overwritten) in init_ocean)
1292       IF ( use_single_reference_value )  THEN
1293          IF (  .NOT.  humidity )  THEN
1294             ref_state(:) = pt_reference
1295          ELSE
1296             ref_state(:) = vpt_reference
1297          ENDIF
1298       ELSE
1299          IF (  .NOT.  humidity )  THEN
1300             ref_state(:) = pt_init(:)
1301          ELSE
1302             ref_state(:) = vpt(:,nys,nxl)
1303          ENDIF
1304       ENDIF
1305
1306!
1307!--    For the moment, vertical velocity is zero
1308       w = 0.0_wp
1309
1310!
1311!--    Initialize array sums (must be defined in first call of pres)
1312       sums = 0.0_wp
1313
1314!
1315!--    In case of iterative solvers, p must get an initial value
1316       IF ( psolver(1:9) == 'multigrid'  .OR.  psolver == 'sor' )  p = 0.0_wp
1317
1318!
1319!--    Treating cloud physics, liquid water content and precipitation amount
1320!--    are zero at beginning of the simulation
1321       IF ( cloud_physics )  THEN
1322          ql = 0.0_wp
1323          qc = 0.0_wp
1324
1325          precipitation_amount = 0.0_wp
1326       ENDIF
1327!
1328!--    Impose vortex with vertical axis on the initial velocity profile
1329       IF ( INDEX( initializing_actions, 'initialize_vortex' ) /= 0 )  THEN
1330          CALL init_rankine
1331       ENDIF
1332
1333!
1334!--    Impose temperature anomaly (advection test only)
1335       IF ( INDEX( initializing_actions, 'initialize_ptanom' ) /= 0 )  THEN
1336          CALL init_pt_anomaly
1337       ENDIF
1338
1339!
1340!--    If required, change the surface temperature at the start of the 3D run
1341       IF ( pt_surface_initial_change /= 0.0_wp )  THEN
1342          pt(nzb,:,:) = pt(nzb,:,:) + pt_surface_initial_change
1343       ENDIF
1344
1345!
1346!--    If required, change the surface humidity/scalar at the start of the 3D
1347!--    run
1348       IF ( humidity  .AND.  q_surface_initial_change /= 0.0_wp )              &
1349          q(nzb,:,:) = q(nzb,:,:) + q_surface_initial_change
1350         
1351       IF ( passive_scalar .AND.  s_surface_initial_change /= 0.0_wp )         &
1352          s(nzb,:,:) = s(nzb,:,:) + s_surface_initial_change
1353       
1354
1355!
1356!--    Initialize old and new time levels.
1357       te_m = 0.0_wp; tpt_m = 0.0_wp; tu_m = 0.0_wp; tv_m = 0.0_wp; tw_m = 0.0_wp
1358       e_p = e; pt_p = pt; u_p = u; v_p = v; w_p = w
1359
1360       IF ( humidity  )  THEN
1361          tq_m = 0.0_wp
1362          q_p = q
1363          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
1364             tqr_m = 0.0_wp
1365             qr_p  = qr
1366             tnr_m = 0.0_wp
1367             nr_p  = nr
1368          ENDIF
1369       ENDIF
1370       
1371       IF ( passive_scalar )  THEN
1372          ts_m = 0.0_wp
1373          s_p  = s
1374       ENDIF       
1375
1376       IF ( ocean )  THEN
1377          tsa_m = 0.0_wp
1378          sa_p  = sa
1379       ENDIF
1380       
1381       CALL location_message( 'finished', .TRUE. )
1382
1383    ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data'  .OR.         &
1384         TRIM( initializing_actions ) == 'cyclic_fill' )                       &
1385    THEN
1386
1387       CALL location_message( 'initializing in case of restart / cyclic_fill', &
1388                              .FALSE. )
1389!
1390!--    When reading data for cyclic fill of 3D prerun data files, read
1391!--    some of the global variables from the restart file which are required
1392!--    for initializing the inflow
1393       IF ( TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
1394
1395          DO  i = 0, io_blocks-1
1396             IF ( i == io_group )  THEN
1397                CALL read_parts_of_var_list
1398                CALL close_file( 13 )
1399             ENDIF
1400#if defined( __parallel )
1401             CALL MPI_BARRIER( comm2d, ierr )
1402#endif
1403          ENDDO
1404
1405       ENDIF
1406
1407!
1408!--    Read binary data from restart file
1409       DO  i = 0, io_blocks-1
1410          IF ( i == io_group )  THEN
1411             CALL read_3d_binary
1412          ENDIF
1413#if defined( __parallel )
1414          CALL MPI_BARRIER( comm2d, ierr )
1415#endif
1416       ENDDO
1417
1418!
1419!--    Initialization of the turbulence recycling method
1420       IF ( TRIM( initializing_actions ) == 'cyclic_fill'  .AND.               &
1421            turbulent_inflow )  THEN
1422!
1423!--       First store the profiles to be used at the inflow.
1424!--       These profiles are the (temporally) and horizontally averaged vertical
1425!--       profiles from the prerun. Alternatively, prescribed profiles
1426!--       for u,v-components can be used.
1427          ALLOCATE( mean_inflow_profiles(nzb:nzt+1,7) )
1428
1429          IF ( use_prescribed_profile_data )  THEN
1430             mean_inflow_profiles(:,1) = u_init            ! u
1431             mean_inflow_profiles(:,2) = v_init            ! v
1432          ELSE
1433             mean_inflow_profiles(:,1) = hom_sum(:,1,0)    ! u
1434             mean_inflow_profiles(:,2) = hom_sum(:,2,0)    ! v
1435          ENDIF
1436          mean_inflow_profiles(:,4) = hom_sum(:,4,0)       ! pt
1437          mean_inflow_profiles(:,5) = hom_sum(:,8,0)       ! e
1438          IF ( humidity )                                                      &
1439             mean_inflow_profiles(:,6) = hom_sum(:,41,0)   ! q
1440          IF ( passive_scalar )                                                &
1441             mean_inflow_profiles(:,7) = hom_sum(:,115,0)   ! s
1442
1443!
1444!--       If necessary, adjust the horizontal flow field to the prescribed
1445!--       profiles
1446          IF ( use_prescribed_profile_data )  THEN
1447             DO  i = nxlg, nxrg
1448                DO  j = nysg, nyng
1449                   DO  k = nzb, nzt+1
1450                      u(k,j,i) = u(k,j,i) - hom_sum(k,1,0) + u_init(k)
1451                      v(k,j,i) = v(k,j,i) - hom_sum(k,2,0) + v_init(k)
1452                   ENDDO
1453                ENDDO
1454             ENDDO
1455          ENDIF
1456
1457!
1458!--       Use these mean profiles at the inflow (provided that Dirichlet
1459!--       conditions are used)
1460          IF ( inflow_l )  THEN
1461             DO  j = nysg, nyng
1462                DO  k = nzb, nzt+1
1463                   u(k,j,nxlg:-1)  = mean_inflow_profiles(k,1)
1464                   v(k,j,nxlg:-1)  = mean_inflow_profiles(k,2)
1465                   w(k,j,nxlg:-1)  = 0.0_wp
1466                   pt(k,j,nxlg:-1) = mean_inflow_profiles(k,4)
1467                   e(k,j,nxlg:-1)  = mean_inflow_profiles(k,5)
1468                   IF ( humidity )                                             &
1469                      q(k,j,nxlg:-1)  = mean_inflow_profiles(k,6)
1470                   IF ( passive_scalar )                                       &
1471                      s(k,j,nxlg:-1)  = mean_inflow_profiles(k,7)                     
1472                ENDDO
1473             ENDDO
1474          ENDIF
1475
1476!
1477!--       Calculate the damping factors to be used at the inflow. For a
1478!--       turbulent inflow the turbulent fluctuations have to be limited
1479!--       vertically because otherwise the turbulent inflow layer will grow
1480!--       in time.
1481          IF ( inflow_damping_height == 9999999.9_wp )  THEN
1482!
1483!--          Default: use the inversion height calculated by the prerun; if
1484!--          this is zero, inflow_damping_height must be explicitly
1485!--          specified.
1486             IF ( hom_sum(nzb+6,pr_palm,0) /= 0.0_wp )  THEN
1487                inflow_damping_height = hom_sum(nzb+6,pr_palm,0)
1488             ELSE
1489                WRITE( message_string, * ) 'inflow_damping_height must be ',   &
1490                     'explicitly specified because&the inversion height ',     &
1491                     'calculated by the prerun is zero.'
1492                CALL message( 'init_3d_model', 'PA0318', 1, 2, 0, 6, 0 )
1493             ENDIF
1494
1495          ENDIF
1496
1497          IF ( inflow_damping_width == 9999999.9_wp )  THEN
1498!
1499!--          Default for the transition range: one tenth of the undamped
1500!--          layer
1501             inflow_damping_width = 0.1_wp * inflow_damping_height
1502
1503          ENDIF
1504
1505          ALLOCATE( inflow_damping_factor(nzb:nzt+1) )
1506
1507          DO  k = nzb, nzt+1
1508
1509             IF ( zu(k) <= inflow_damping_height )  THEN
1510                inflow_damping_factor(k) = 1.0_wp
1511             ELSEIF ( zu(k) <= ( inflow_damping_height + inflow_damping_width ) )  THEN
1512                inflow_damping_factor(k) = 1.0_wp -                            &
1513                                           ( zu(k) - inflow_damping_height ) / &
1514                                           inflow_damping_width
1515             ELSE
1516                inflow_damping_factor(k) = 0.0_wp
1517             ENDIF
1518
1519          ENDDO
1520
1521       ENDIF
1522
1523!
1524!--    Inside buildings set velocities and TKE back to zero
1525       IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND.                &
1526            topography /= 'flat' )  THEN
1527!
1528!--       Inside buildings set velocities and TKE back to zero.
1529!--       Other scalars (pt, q, s, km, kh, p, sa, ...) are ignored at present,
1530!--       maybe revise later.
1531          DO  i = nxlg, nxrg
1532             DO  j = nysg, nyng
1533                u  (nzb:nzb_u_inner(j,i),j,i)   = 0.0_wp
1534                v  (nzb:nzb_v_inner(j,i),j,i)   = 0.0_wp
1535                w  (nzb:nzb_w_inner(j,i),j,i)   = 0.0_wp
1536                e  (nzb:nzb_w_inner(j,i),j,i)   = 0.0_wp
1537                tu_m(nzb:nzb_u_inner(j,i),j,i)  = 0.0_wp
1538                tv_m(nzb:nzb_v_inner(j,i),j,i)  = 0.0_wp
1539                tw_m(nzb:nzb_w_inner(j,i),j,i)  = 0.0_wp
1540                te_m(nzb:nzb_w_inner(j,i),j,i)  = 0.0_wp
1541                tpt_m(nzb:nzb_w_inner(j,i),j,i) = 0.0_wp
1542             ENDDO
1543          ENDDO
1544
1545       ENDIF
1546
1547!
1548!--    Calculate initial temperature field and other constants used in case
1549!--    of a sloping surface
1550       IF ( sloping_surface )  CALL init_slope
1551
1552!
1553!--    Initialize new time levels (only done in order to set boundary values
1554!--    including ghost points)
1555       e_p = e; pt_p = pt; u_p = u; v_p = v; w_p = w
1556       IF ( humidity )  THEN
1557          q_p = q
1558          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
1559             qr_p = qr
1560             nr_p = nr
1561          ENDIF
1562       ENDIF
1563       IF ( passive_scalar )  s_p  = s
1564       IF ( ocean          )  sa_p = sa
1565
1566!
1567!--    Allthough tendency arrays are set in prognostic_equations, they have
1568!--    have to be predefined here because they are used (but multiplied with 0)
1569!--    there before they are set.
1570       te_m = 0.0_wp; tpt_m = 0.0_wp; tu_m = 0.0_wp; tv_m = 0.0_wp; tw_m = 0.0_wp
1571       IF ( humidity )  THEN
1572          tq_m = 0.0_wp
1573          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
1574             tqr_m = 0.0_wp
1575             tnr_m = 0.0_wp
1576          ENDIF
1577       ENDIF
1578       IF ( passive_scalar )  ts_m  = 0.0_wp
1579       IF ( ocean          )  tsa_m = 0.0_wp
1580
1581       CALL location_message( 'finished', .TRUE. )
1582
1583    ELSE
1584!
1585!--    Actually this part of the programm should not be reached
1586       message_string = 'unknown initializing problem'
1587       CALL message( 'init_3d_model', 'PA0193', 1, 2, 0, 6, 0 )
1588    ENDIF
1589
1590
1591    IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1592!
1593!--    Initialize old timelevels needed for radiation boundary conditions
1594       IF ( outflow_l )  THEN
1595          u_m_l(:,:,:) = u(:,:,1:2)
1596          v_m_l(:,:,:) = v(:,:,0:1)
1597          w_m_l(:,:,:) = w(:,:,0:1)
1598       ENDIF
1599       IF ( outflow_r )  THEN
1600          u_m_r(:,:,:) = u(:,:,nx-1:nx)
1601          v_m_r(:,:,:) = v(:,:,nx-1:nx)
1602          w_m_r(:,:,:) = w(:,:,nx-1:nx)
1603       ENDIF
1604       IF ( outflow_s )  THEN
1605          u_m_s(:,:,:) = u(:,0:1,:)
1606          v_m_s(:,:,:) = v(:,1:2,:)
1607          w_m_s(:,:,:) = w(:,0:1,:)
1608       ENDIF
1609       IF ( outflow_n )  THEN
1610          u_m_n(:,:,:) = u(:,ny-1:ny,:)
1611          v_m_n(:,:,:) = v(:,ny-1:ny,:)
1612          w_m_n(:,:,:) = w(:,ny-1:ny,:)
1613       ENDIF
1614       
1615    ENDIF
1616
1617!
1618!-- Calculate the initial volume flow at the right and north boundary
1619    IF ( conserve_volume_flow )  THEN
1620
1621       IF ( use_prescribed_profile_data )  THEN
1622
1623          volume_flow_initial_l = 0.0_wp
1624          volume_flow_area_l    = 0.0_wp
1625
1626          IF ( nxr == nx )  THEN
1627             DO  j = nys, nyn
1628                DO  k = nzb_u_inner(j,nx)+1, nzt
1629                   volume_flow_initial_l(1) = volume_flow_initial_l(1) +       &
1630                                              u_init(k) * dzw(k)
1631                   volume_flow_area_l(1)    = volume_flow_area_l(1) + dzw(k)
1632                ENDDO
1633             ENDDO
1634          ENDIF
1635         
1636          IF ( nyn == ny )  THEN
1637             DO  i = nxl, nxr
1638                DO  k = nzb_v_inner(ny,i)+1, nzt
1639                   volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
1640                                              v_init(k) * dzw(k)
1641                   volume_flow_area_l(2)    = volume_flow_area_l(2) + dzw(k)
1642                ENDDO
1643             ENDDO
1644          ENDIF
1645
1646#if defined( __parallel )
1647          CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&
1648                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
1649          CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1),      &
1650                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
1651
1652#else
1653          volume_flow_initial = volume_flow_initial_l
1654          volume_flow_area    = volume_flow_area_l
1655#endif 
1656
1657       ELSEIF ( TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
1658
1659          volume_flow_initial_l = 0.0_wp
1660          volume_flow_area_l    = 0.0_wp
1661
1662          IF ( nxr == nx )  THEN
1663             DO  j = nys, nyn
1664                DO  k = nzb_u_inner(j,nx)+1, nzt
1665                   volume_flow_initial_l(1) = volume_flow_initial_l(1) +       &
1666                                              hom_sum(k,1,0) * dzw(k)
1667                   volume_flow_area_l(1)    = volume_flow_area_l(1) + dzw(k)
1668                ENDDO
1669             ENDDO
1670          ENDIF
1671         
1672          IF ( nyn == ny )  THEN
1673             DO  i = nxl, nxr
1674                DO  k = nzb_v_inner(ny,i)+1, nzt
1675                   volume_flow_initial_l(2) = volume_flow_initial_l(2) +       &
1676                                              hom_sum(k,2,0) * dzw(k)
1677                   volume_flow_area_l(2)    = volume_flow_area_l(2) + dzw(k)
1678                ENDDO
1679             ENDDO
1680          ENDIF
1681
1682#if defined( __parallel )
1683          CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&
1684                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
1685          CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1),      &
1686                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
1687
1688#else
1689          volume_flow_initial = volume_flow_initial_l
1690          volume_flow_area    = volume_flow_area_l
1691#endif 
1692
1693       ELSEIF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1694
1695          volume_flow_initial_l = 0.0_wp
1696          volume_flow_area_l    = 0.0_wp
1697
1698          IF ( nxr == nx )  THEN
1699             DO  j = nys, nyn
1700                DO  k = nzb_u_inner(j,nx)+1, nzt
1701                   volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
1702                                              u(k,j,nx) * dzw(k)
1703                   volume_flow_area_l(1)    = volume_flow_area_l(1) + dzw(k)
1704                ENDDO
1705             ENDDO
1706          ENDIF
1707         
1708          IF ( nyn == ny )  THEN
1709             DO  i = nxl, nxr
1710                DO  k = nzb_v_inner(ny,i)+1, nzt
1711                   volume_flow_initial_l(2) = volume_flow_initial_l(2) +       &
1712                                              v(k,ny,i) * dzw(k)
1713                   volume_flow_area_l(2)    = volume_flow_area_l(2) + dzw(k)
1714                ENDDO
1715             ENDDO
1716          ENDIF
1717
1718#if defined( __parallel )
1719          CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&
1720                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
1721          CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1),      &
1722                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
1723
1724#else
1725          volume_flow_initial = volume_flow_initial_l
1726          volume_flow_area    = volume_flow_area_l
1727#endif 
1728
1729       ENDIF
1730
1731!
1732!--    In case of 'bulk_velocity' mode, volume_flow_initial is calculated
1733!--    from u|v_bulk instead
1734       IF ( TRIM( conserve_volume_flow_mode ) == 'bulk_velocity' )  THEN
1735          volume_flow_initial(1) = u_bulk * volume_flow_area(1)
1736          volume_flow_initial(2) = v_bulk * volume_flow_area(2)
1737       ENDIF
1738
1739    ENDIF
1740
1741!
1742!-- Initialize quantities for special advections schemes
1743    CALL init_advec
1744
1745!
1746!-- Impose random perturbation on the horizontal velocity field and then
1747!-- remove the divergences from the velocity field at the initial stage
1748    IF ( create_disturbances  .AND.  disturbance_energy_limit /= 0.0_wp  .AND. &
1749         TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
1750         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
1751
1752       CALL location_message( 'creating initial disturbances', .FALSE. )
1753       CALL disturb_field( nzb_u_inner, tend, u )
1754       CALL disturb_field( nzb_v_inner, tend, v )
1755       CALL location_message( 'finished', .TRUE. )
1756
1757       CALL location_message( 'calling pressure solver', .FALSE. )
1758       n_sor = nsor_ini
1759       !$acc data copyin( d, ddzu, ddzw, nzb_s_inner, nzb_u_inner )            &
1760       !$acc      copyin( nzb_v_inner, nzb_w_inner, p, rflags_s_inner, tend )  &
1761       !$acc      copyin( weight_pres, weight_substep )                        &
1762       !$acc      copy( tri, tric, u, v, w )
1763       CALL pres
1764       !$acc end data
1765       n_sor = nsor
1766       CALL location_message( 'finished', .TRUE. )
1767
1768    ENDIF
1769
1770!
1771!-- If required, initialize quantities needed for the plant canopy model
1772    IF ( plant_canopy )  THEN
1773       CALL location_message( 'initializing plant canopy model', .FALSE. )   
1774       CALL pcm_init
1775       CALL location_message( 'finished', .TRUE. )
1776    ENDIF
1777
1778!
1779!-- If required, initialize dvrp-software
1780    IF ( dt_dvrp /= 9999999.9_wp )  CALL init_dvrp
1781
1782    IF ( ocean )  THEN
1783!
1784!--    Initialize quantities needed for the ocean model
1785       CALL init_ocean
1786
1787    ELSE
1788!
1789!--    Initialize quantities for handling cloud physics
1790!--    This routine must be called before lpm_init, because
1791!--    otherwise, array pt_d_t, needed in data_output_dvrp (called by
1792!--    lpm_init) is not defined.
1793       CALL init_cloud_physics
1794!
1795!--    Initialize bulk cloud microphysics
1796       CALL microphysics_init
1797    ENDIF
1798
1799!
1800!-- If required, initialize particles
1801    IF ( particle_advection )  CALL lpm_init
1802
1803!
1804!-- If required, initialize quantities needed for the LSM
1805    IF ( land_surface )  THEN
1806       CALL location_message( 'initializing land surface model', .FALSE. )
1807       CALL lsm_init
1808       CALL location_message( 'finished', .TRUE. )
1809    ENDIF
1810
1811!
1812!-- Initialize surface layer (done after LSM as roughness length are required
1813!-- for initialization
1814    IF ( constant_flux_layer )  THEN
1815       CALL location_message( 'initializing surface layer', .FALSE. )
1816       CALL init_surface_layer_fluxes
1817       CALL location_message( 'finished', .TRUE. )
1818    ENDIF
1819
1820!
1821!-- If required, initialize radiation model
1822    IF ( radiation )  THEN
1823       CALL location_message( 'initializing radiation model', .FALSE. )
1824       CALL radiation_init
1825       CALL location_message( 'finished', .TRUE. )
1826    ENDIF
1827
1828!
1829!-- If required, initialize urban surface model
1830    IF ( urban_surface )  THEN
1831       CALL location_message( 'initializing urban surface model', .FALSE. )
1832       CALL usm_init_urban_surface
1833       CALL location_message( 'finished', .TRUE. )
1834    ENDIF
1835
1836!
1837!-- If required, initialize quantities needed for the wind turbine model
1838    IF ( wind_turbine )  THEN
1839       CALL location_message( 'initializing wind turbine model', .FALSE. )
1840       CALL wtm_init
1841       CALL location_message( 'finished', .TRUE. )
1842    ENDIF
1843
1844
1845!
1846!-- Initialize the ws-scheme.   
1847    IF ( ws_scheme_sca .OR. ws_scheme_mom )  CALL ws_init       
1848
1849!
1850!-- Setting weighting factors for calculation of perturbation pressure
1851!-- and turbulent quantities from the RK substeps
1852    IF ( TRIM(timestep_scheme) == 'runge-kutta-3' )  THEN      ! for RK3-method
1853
1854       weight_substep(1) = 1._wp/6._wp
1855       weight_substep(2) = 3._wp/10._wp
1856       weight_substep(3) = 8._wp/15._wp
1857
1858       weight_pres(1)    = 1._wp/3._wp
1859       weight_pres(2)    = 5._wp/12._wp
1860       weight_pres(3)    = 1._wp/4._wp
1861
1862    ELSEIF ( TRIM(timestep_scheme) == 'runge-kutta-2' )  THEN  ! for RK2-method
1863
1864       weight_substep(1) = 1._wp/2._wp
1865       weight_substep(2) = 1._wp/2._wp
1866         
1867       weight_pres(1)    = 1._wp/2._wp
1868       weight_pres(2)    = 1._wp/2._wp       
1869
1870    ELSE                                     ! for Euler-method
1871
1872       weight_substep(1) = 1.0_wp     
1873       weight_pres(1)    = 1.0_wp                   
1874
1875    ENDIF
1876
1877!
1878!-- Initialize Rayleigh damping factors
1879    rdf    = 0.0_wp
1880    rdf_sc = 0.0_wp
1881    IF ( rayleigh_damping_factor /= 0.0_wp )  THEN
1882       IF (  .NOT.  ocean )  THEN
1883          DO  k = nzb+1, nzt
1884             IF ( zu(k) >= rayleigh_damping_height )  THEN
1885                rdf(k) = rayleigh_damping_factor *                             &
1886                      ( SIN( pi * 0.5_wp * ( zu(k) - rayleigh_damping_height ) &
1887                             / ( zu(nzt) - rayleigh_damping_height ) )         &
1888                      )**2
1889             ENDIF
1890          ENDDO
1891       ELSE
1892          DO  k = nzt, nzb+1, -1
1893             IF ( zu(k) <= rayleigh_damping_height )  THEN
1894                rdf(k) = rayleigh_damping_factor *                             &
1895                      ( SIN( pi * 0.5_wp * ( rayleigh_damping_height - zu(k) ) &
1896                             / ( rayleigh_damping_height - zu(nzb+1) ) )       &
1897                      )**2
1898             ENDIF
1899          ENDDO
1900       ENDIF
1901    ENDIF
1902    IF ( scalar_rayleigh_damping )  rdf_sc = rdf
1903
1904!
1905!-- Initialize the starting level and the vertical smoothing factor used for
1906!-- the external pressure gradient
1907    dp_smooth_factor = 1.0_wp
1908    IF ( dp_external )  THEN
1909!
1910!--    Set the starting level dp_level_ind_b only if it has not been set before
1911!--    (e.g. in init_grid).
1912       IF ( dp_level_ind_b == 0 )  THEN
1913          ind_array = MINLOC( ABS( dp_level_b - zu ) )
1914          dp_level_ind_b = ind_array(1) - 1 + nzb
1915                                        ! MINLOC uses lower array bound 1
1916       ENDIF
1917       IF ( dp_smooth )  THEN
1918          dp_smooth_factor(:dp_level_ind_b) = 0.0_wp
1919          DO  k = dp_level_ind_b+1, nzt
1920             dp_smooth_factor(k) = 0.5_wp * ( 1.0_wp + SIN( pi *               &
1921                        ( REAL( k - dp_level_ind_b, KIND=wp ) /                &
1922                          REAL( nzt - dp_level_ind_b, KIND=wp ) - 0.5_wp ) ) )
1923          ENDDO
1924       ENDIF
1925    ENDIF
1926
1927!
1928!-- Initialize damping zone for the potential temperature in case of
1929!-- non-cyclic lateral boundaries. The damping zone has the maximum value
1930!-- at the inflow boundary and decreases to zero at pt_damping_width.
1931    ptdf_x = 0.0_wp
1932    ptdf_y = 0.0_wp
1933    IF ( bc_lr_dirrad )  THEN
1934       DO  i = nxl, nxr
1935          IF ( ( i * dx ) < pt_damping_width )  THEN
1936             ptdf_x(i) = pt_damping_factor * ( SIN( pi * 0.5_wp *              &
1937                            REAL( pt_damping_width - i * dx, KIND=wp ) / (     &
1938                            REAL( pt_damping_width, KIND=wp ) ) ) )**2 
1939          ENDIF
1940       ENDDO
1941    ELSEIF ( bc_lr_raddir )  THEN
1942       DO  i = nxl, nxr
1943          IF ( ( i * dx ) > ( nx * dx - pt_damping_width ) )  THEN
1944             ptdf_x(i) = pt_damping_factor *                                   &
1945                         SIN( pi * 0.5_wp *                                    &
1946                                 ( ( i - nx ) * dx + pt_damping_width ) /      &
1947                                 REAL( pt_damping_width, KIND=wp ) )**2
1948          ENDIF
1949       ENDDO
1950    ELSEIF ( bc_ns_dirrad )  THEN
1951       DO  j = nys, nyn
1952          IF ( ( j * dy ) > ( ny * dy - pt_damping_width ) )  THEN
1953             ptdf_y(j) = pt_damping_factor *                                   &
1954                         SIN( pi * 0.5_wp *                                    &
1955                                 ( ( j - ny ) * dy + pt_damping_width ) /      &
1956                                 REAL( pt_damping_width, KIND=wp ) )**2
1957          ENDIF
1958       ENDDO
1959    ELSEIF ( bc_ns_raddir )  THEN
1960       DO  j = nys, nyn
1961          IF ( ( j * dy ) < pt_damping_width )  THEN
1962             ptdf_y(j) = pt_damping_factor *                                   &
1963                         SIN( pi * 0.5_wp *                                    &
1964                                ( pt_damping_width - j * dy ) /                &
1965                                REAL( pt_damping_width, KIND=wp ) )**2
1966          ENDIF
1967       ENDDO
1968    ENDIF
1969
1970!
1971!-- Pre-set masks for regional statistics. Default is the total model domain.
1972!-- Ghost points are excluded because counting values at the ghost boundaries
1973!-- would bias the statistics
1974    rmask = 1.0_wp
1975    rmask(:,nxlg:nxl-1,:) = 0.0_wp;  rmask(:,nxr+1:nxrg,:) = 0.0_wp
1976    rmask(nysg:nys-1,:,:) = 0.0_wp;  rmask(nyn+1:nyng,:,:) = 0.0_wp
1977
1978!
1979!-- User-defined initializing actions. Check afterwards, if maximum number
1980!-- of allowed timeseries is exceeded
1981    CALL user_init
1982
1983    IF ( dots_num > dots_max )  THEN
1984       WRITE( message_string, * ) 'number of time series quantities exceeds',  &
1985                                  ' its maximum of dots_max = ', dots_max,     &
1986                                  ' &Please increase dots_max in modules.f90.'
1987       CALL message( 'init_3d_model', 'PA0194', 1, 2, 0, 6, 0 )   
1988    ENDIF
1989
1990!
1991!-- Input binary data file is not needed anymore. This line must be placed
1992!-- after call of user_init!
1993    CALL close_file( 13 )
1994
1995!
1996!-- Compute total sum of active mask grid points
1997!-- and the mean surface level height for each statistic region
1998!-- ngp_2dh: number of grid points of a horizontal cross section through the
1999!--          total domain
2000!-- ngp_3d:  number of grid points of the total domain
2001    ngp_2dh_outer_l   = 0
2002    ngp_2dh_outer     = 0
2003    ngp_2dh_s_inner_l = 0
2004    ngp_2dh_s_inner   = 0
2005    ngp_2dh_l         = 0
2006    ngp_2dh           = 0
2007    ngp_3d_inner_l    = 0.0_wp
2008    ngp_3d_inner      = 0
2009    ngp_3d            = 0
2010    ngp_sums          = ( nz + 2 ) * ( pr_palm + max_pr_user )
2011
2012    mean_surface_level_height   = 0.0_wp
2013    mean_surface_level_height_l = 0.0_wp
2014
2015    DO  sr = 0, statistic_regions
2016       DO  i = nxl, nxr
2017          DO  j = nys, nyn
2018             IF ( rmask(j,i,sr) == 1.0_wp )  THEN
2019!
2020!--             All xy-grid points
2021                ngp_2dh_l(sr) = ngp_2dh_l(sr) + 1
2022                mean_surface_level_height_l(sr) = mean_surface_level_height_l(sr) &
2023                                                  + zw(nzb_s_inner(j,i))
2024!
2025!--             xy-grid points above topography
2026                DO  k = nzb_s_outer(j,i), nz + 1
2027                   ngp_2dh_outer_l(k,sr) = ngp_2dh_outer_l(k,sr) + 1
2028                ENDDO
2029                DO  k = nzb_s_inner(j,i), nz + 1
2030                   ngp_2dh_s_inner_l(k,sr) = ngp_2dh_s_inner_l(k,sr) + 1
2031                ENDDO
2032!
2033!--             All grid points of the total domain above topography
2034                ngp_3d_inner_l(sr) = ngp_3d_inner_l(sr)                        &
2035                                     + ( nz - nzb_s_inner(j,i) + 2 )
2036             ENDIF
2037          ENDDO
2038       ENDDO
2039    ENDDO
2040
2041    sr = statistic_regions + 1
2042#if defined( __parallel )
2043    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2044    CALL MPI_ALLREDUCE( ngp_2dh_l(0), ngp_2dh(0), sr, MPI_INTEGER, MPI_SUM,    &
2045                        comm2d, ierr )
2046    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2047    CALL MPI_ALLREDUCE( ngp_2dh_outer_l(0,0), ngp_2dh_outer(0,0), (nz+2)*sr,   &
2048                        MPI_INTEGER, MPI_SUM, comm2d, ierr )
2049    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2050    CALL MPI_ALLREDUCE( ngp_2dh_s_inner_l(0,0), ngp_2dh_s_inner(0,0),          &
2051                        (nz+2)*sr, MPI_INTEGER, MPI_SUM, comm2d, ierr )
2052    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2053    CALL MPI_ALLREDUCE( ngp_3d_inner_l(0), ngp_3d_inner_tmp(0), sr, MPI_REAL,  &
2054                        MPI_SUM, comm2d, ierr )
2055    ngp_3d_inner = INT( ngp_3d_inner_tmp, KIND = SELECTED_INT_KIND( 18 ) )
2056    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2057    CALL MPI_ALLREDUCE( mean_surface_level_height_l(0),                        &
2058                        mean_surface_level_height(0), sr, MPI_REAL,            &
2059                        MPI_SUM, comm2d, ierr )
2060    mean_surface_level_height = mean_surface_level_height / REAL( ngp_2dh )
2061#else
2062    ngp_2dh         = ngp_2dh_l
2063    ngp_2dh_outer   = ngp_2dh_outer_l
2064    ngp_2dh_s_inner = ngp_2dh_s_inner_l
2065    ngp_3d_inner    = INT( ngp_3d_inner_l, KIND = SELECTED_INT_KIND( 18 ) )
2066    mean_surface_level_height = mean_surface_level_height_l / REAL( ngp_2dh_l )
2067#endif
2068
2069    ngp_3d = INT ( ngp_2dh, KIND = SELECTED_INT_KIND( 18 ) ) * &
2070             INT ( (nz + 2 ), KIND = SELECTED_INT_KIND( 18 ) )
2071
2072!
2073!-- Set a lower limit of 1 in order to avoid zero divisions in flow_statistics,
2074!-- buoyancy, etc. A zero value will occur for cases where all grid points of
2075!-- the respective subdomain lie below the surface topography
2076    ngp_2dh_outer   = MAX( 1, ngp_2dh_outer(:,:)   ) 
2077    ngp_3d_inner    = MAX( INT(1, KIND = SELECTED_INT_KIND( 18 )),             &
2078                           ngp_3d_inner(:) )
2079    ngp_2dh_s_inner = MAX( 1, ngp_2dh_s_inner(:,:) ) 
2080
2081    DEALLOCATE( mean_surface_level_height_l, ngp_2dh_l, ngp_2dh_outer_l,       &
2082                ngp_3d_inner_l, ngp_3d_inner_tmp )
2083
2084    CALL location_message( 'leaving init_3d_model', .TRUE. )
2085
2086 END SUBROUTINE init_3d_model
Note: See TracBrowser for help on using the repository browser.