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

Last change on this file since 4131 was 4131, checked in by monakurppa, 2 years ago

Several changes in the salsa aerosol module:

  • Add "salsa_" before each salsa output variable
  • Add a possibility to output the number (salsa_N_UFP) and mass concentration (salsa_PM0.1) of ultrafine particles, i.e. particles with a diameter smaller than 100 nm
  • Implement aerosol emission mode "parameterized" which is based on the street type (similar to the chemistry module).
  • Remove unnecessary nucleation subroutines.
  • Add the z-dimension for gaseous emissions to correspond the implementation in the chemistry module
  • Property svn:keywords set to Id
  • Property svn:mergeinfo set to False
    /palm/branches/chemistry/SOURCE/init_3d_model.f902047-3190,​3218-3297
    /palm/branches/forwind/SOURCE/init_3d_model.f901564-1913
    /palm/branches/mosaik_M2/init_3d_model.f902360-3471
    /palm/branches/palm4u/SOURCE/init_3d_model.f902540-2692
    /palm/branches/rans/SOURCE/init_3d_model.f902078-3128
    /palm/branches/resler/SOURCE/init_3d_model.f902023-3605
    /palm/branches/salsa/SOURCE/init_3d_model.f902503-3581
File size: 90.6 KB
Line 
1!> @file init_3d_model.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
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-2019 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: init_3d_model.f90 4131 2019-08-02 11:06:18Z monakurppa $
27! Allocate sums and sums_l to allow profile output for salsa variables.
28!
29! 4130 2019-08-01 13:04:13Z suehring
30! Effectively reduce 3D initialization to 1D initial profiles. This is because
31! 3D initialization produces structures in the w-component that are correlated
32! with the processor grid for some unknown reason 
33!
34! 4090 2019-07-11 15:06:47Z Giersch
35! Unused variables removed
36!
37! 4088 2019-07-11 13:57:56Z Giersch
38! Pressure and density profile calculations revised using basic functions
39!
40! 4048 2019-06-21 21:00:21Z knoop
41! Further modularization of particle code components
42!
43! 4017 2019-06-06 12:16:46Z schwenkel
44! Convert most location messages to debug messages to reduce output in
45! job logfile to a minimum
46!
47!
48! unused variable removed
49!
50! 3937 2019-04-29 15:09:07Z suehring
51! Move initialization of synthetic turbulence generator behind initialization
52! of offline nesting. Remove call for stg_adjust, as this is now already done
53! in stg_init.
54!
55! 3900 2019-04-16 15:17:43Z suehring
56! Fix problem with LOD = 2 initialization
57!
58! 3885 2019-04-11 11:29:34Z kanani
59! Changes related to global restructuring of location messages and introduction
60! of additional debug messages
61!
62! 3849 2019-04-01 16:35:16Z knoop
63! Move initialization of rmask before initializing user_init_arrays
64!
65! 3711 2019-01-31 13:44:26Z knoop
66! Introduced module_interface_init_checks for post-init checks in modules
67!
68! 3700 2019-01-26 17:03:42Z knoop
69! Some interface calls moved to module_interface + cleanup
70!
71! 3648 2019-01-02 16:35:46Z suehring
72! Rename subroutines for surface-data output
73!
74! 3636 2018-12-19 13:48:34Z raasch
75! nopointer option removed
76!
77! 3609 2018-12-07 13:37:59Z suehring
78! Furhter correction in initialization of surfaces in cyclic-fill case
79!
80! 3608 2018-12-07 12:59:57Z suehring
81! Bugfix in initialization of surfaces in cyclic-fill case
82!
83! 3589 2018-11-30 15:09:51Z suehring
84! Move the control parameter "salsa" from salsa_mod to control_parameters
85! (M. Kurppa)
86!
87! 3582 2018-11-29 19:16:36Z suehring
88! Bugfix in initialization of turbulence generator
89!
90! 3569 2018-11-27 17:03:40Z kanani
91! dom_dwd_user, Schrempf:
92! Remove uv exposure model code, this is now part of biometeorology_mod,
93! remove bio_init_arrays.
94!
95! 3547 2018-11-21 13:21:24Z suehring
96! variables documented
97!
98! 3525 2018-11-14 16:06:14Z kanani
99! Changes related to clean-up of biometeorology (dom_dwd_user)
100!
101! 3524 2018-11-14 13:36:44Z raasch
102! preprocessor directive added to avoid the compiler to complain about unused
103! variable
104!
105! 3473 2018-10-30 20:50:15Z suehring
106! Add virtual measurement module
107!
108! 3472 2018-10-30 20:43:50Z suehring
109! Add indoor model (kanani, srissman, tlang)
110!
111! 3467 2018-10-30 19:05:21Z suehring
112! Implementation of a new aerosol module salsa.
113!
114! 3458 2018-10-30 14:51:23Z kanani
115! from chemistry branch r3443, basit:
116! bug fixed in sums and sums_l for chemistry profile output
117!
118! 3448 2018-10-29 18:14:31Z kanani
119! Add biometeorology
120!
121! 3421 2018-10-24 18:39:32Z gronemeier
122! Initialize surface data output
123!
124! 3415 2018-10-24 11:57:50Z suehring
125! Set bottom boundary condition for geostrophic wind components in inifor
126! initialization
127!
128! 3347 2018-10-15 14:21:08Z suehring
129! - Separate offline nesting from large_scale_nudging_mod
130! - Improve the synthetic turbulence generator
131!
132! 3298 2018-10-02 12:21:11Z kanani
133! Minor formatting (kanani)
134! Added CALL to the chem_emissions module (Russo)
135!
136! 3294 2018-10-01 02:37:10Z raasch
137! allocate and set stokes drift velocity profiles
138!
139! 3298 2018-10-02 12:21:11Z kanani
140! Minor formatting (kanani)
141! Added CALL to the chem_emissions module (Russo)
142!
143! 3294 2018-10-01 02:37:10Z raasch
144! changes concerning modularization of ocean option
145!
146! 3289 2018-09-28 10:23:58Z suehring
147! Introduce module parameter for number of inflow profiles
148!
149! 3288 2018-09-28 10:23:08Z suehring
150! Modularization of all bulk cloud physics code components
151!
152! 3241 2018-09-12 15:02:00Z raasch
153! unused variables removed
154!
155! 3234 2018-09-07 13:46:58Z schwenkel
156! The increase of dots_num in case of radiation or land surface model must
157! be done before user_init is called
158!
159! 3183 2018-07-27 14:25:55Z suehring
160! Revise Inifor initialization
161!
162! 3182 2018-07-27 13:36:03Z suehring
163! Added multi agent system
164!
165! 3129 2018-07-16 07:45:13Z gronemeier
166! Move initialization call for nudging and 1D/3D offline nesting.
167! Revise initialization with inifor data.
168!
169! 3045 2018-05-28 07:55:41Z Giersch
170! Error messages revised
171!
172! 3045 2018-05-28 07:55:41Z Giersch
173! Error messages revised
174!
175! 3042 2018-05-25 10:44:37Z schwenkel
176! Changed the name specific humidity to mixing ratio
177!
178! 3040 2018-05-25 10:22:08Z schwenkel
179! Add option to initialize warm air bubble close to surface
180!
181! 3014 2018-05-09 08:42:38Z maronga
182! Bugfix: initialization of ts_value missing
183!
184! 3011 2018-05-07 14:38:42Z schwenkel
185! removed redundant if statement
186!
187! 3004 2018-04-27 12:33:25Z Giersch
188! precipitation_rate removed
189!
190! 2995 2018-04-19 12:13:16Z Giersch
191! CALL radiation_control is not necessary during initialization because
192! calculation of radiative fluxes at model start is done in radiation_init
193! in any case
194!
195! 2977 2018-04-17 10:27:57Z kanani
196! Implement changes from branch radiation (r2948-2971) with minor modifications
197! (moh.hefny):
198! - set radiation_interactions according to the existence of urban/land vertical
199!   surfaces and trees to activiate RTM
200! - set average_radiation to TRUE if RTM is activiated
201!
202! 2938 2018-03-27 15:52:42Z suehring
203! - Revise Inifor initialization for geostrophic wind components
204! - Initialize synthetic turbulence generator in case of Inifor initialization 
205!
206! 2936 2018-03-27 14:49:27Z suehring
207! Synchronize parent and child models after initialization.
208! Remove obsolete masking of topography grid points for Runge-Kutta weighted
209! tendency arrays.
210!
211! 2920 2018-03-22 11:22:01Z kanani
212! Add call for precalculating apparent solar positions (moh.hefny)
213!
214! 2906 2018-03-19 08:56:40Z Giersch
215! The variables read/write_svf_on_init have been removed. Instead ENVIRONMENT
216! variables read/write_svf have been introduced. Location_message has been
217! added.
218!
219! 2894 2018-03-15 09:17:58Z Giersch
220! Renamed routines with respect to reading restart data, file 13 is closed in
221! rrd_read_parts_of_global now
222!
223! 2867 2018-03-09 09:40:23Z suehring
224! Further bugfix concerning call of user_init.
225!
226! 2864 2018-03-08 11:57:45Z suehring
227! Bugfix, move call of user_init in front of initialization of grid-point
228! arrays
229!
230! 2817 2018-02-19 16:32:21Z knoop
231! Preliminary gust module interface implemented
232!
233! 2776 2018-01-31 10:44:42Z Giersch
234! Variable use_synthetic_turbulence_generator has been abbreviated
235!
236! 2766 2018-01-22 17:17:47Z kanani
237! Removed preprocessor directive __chem
238!
239! 2758 2018-01-17 12:55:21Z suehring
240! In case of spinup of land- and urban-surface model, do not mask wind velocity
241! at first computational grid level
242!
243! 2746 2018-01-15 12:06:04Z suehring
244! Move flag plant canopy to modules
245!
246! 2718 2018-01-02 08:49:38Z maronga
247! Corrected "Former revisions" section
248!
249! 2705 2017-12-18 11:26:23Z maronga
250! Bugfix for reading initial profiles from ls/nuding file
251!
252! 2701 2017-12-15 15:40:50Z suehring
253! Changes from last commit documented
254!
255! 2700 2017-12-15 14:12:35Z suehring
256! Bugfix, missing initialization of surface attributes in case of
257! inifor-initialization branch
258!
259! 2698 2017-12-14 18:46:24Z suehring
260! Bugfix in get_topography_top_index
261!
262! 2696 2017-12-14 17:12:51Z kanani
263! Change in file header (GPL part)
264! Implementation of uv exposure model (FK)
265! Moved initialisation of diss, e, kh, km to turbulence_closure_mod (TG)
266! Added chemical emissions (FK)
267! Initialize masking arrays and number-of-grid-points arrays before initialize
268! LSM, USM and radiation module
269! Initialization with inifor (MS)
270!
271! 2618 2017-11-16 15:37:30Z suehring
272! Reorder calls of init_surfaces.
273!
274! 2564 2017-10-19 15:56:56Z Giersch
275! Variable wind_turbine was added to control_parameters.
276!
277! 2550 2017-10-16 17:12:01Z boeske
278! Modifications to cyclic fill method and turbulence recycling method in case of
279! complex terrain simulations
280!
281! 2513 2017-10-04 09:24:39Z kanani
282! Bugfix in storing initial scalar profile (wrong index)
283!
284! 2350 2017-08-15 11:48:26Z kanani
285! Bugfix in nopointer version
286!
287! 2339 2017-08-07 13:55:26Z gronemeier
288! corrected timestamp in header
289!
290! 2338 2017-08-07 12:15:38Z gronemeier
291! Modularize 1D model
292!
293! 2329 2017-08-03 14:24:56Z knoop
294! Removed temporary bugfix (r2327) as bug is properly resolved by this revision
295!
296! 2327 2017-08-02 07:40:57Z maronga
297! Temporary bugfix
298!
299! 2320 2017-07-21 12:47:43Z suehring
300! Modularize large-scale forcing and nudging
301!
302! 2292 2017-06-20 09:51:42Z schwenkel
303! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
304! includes two more prognostic equations for cloud drop concentration (nc) 
305! and cloud water content (qc).
306!
307! 2277 2017-06-12 10:47:51Z kanani
308! Removed unused variable sums_up_fraction_l
309!
310! 2270 2017-06-09 12:18:47Z maronga
311! dots_num must be increased when LSM and/or radiation is used
312!
313! 2259 2017-06-08 09:09:11Z gronemeier
314! Implemented synthetic turbulence generator
315!
316! 2252 2017-06-07 09:35:37Z knoop
317! rho_air now depending on surface_pressure even in Boussinesq mode
318!
319! 2233 2017-05-30 18:08:54Z suehring
320!
321! 2232 2017-05-30 17:47:52Z suehring
322! Adjustments to new topography and surface concept:
323!   - Modify passed parameters for disturb_field
324!   - Topography representation via flags
325!   - Remove unused arrays.
326!   - Move initialization of surface-related quantities to surface_mod
327!
328! 2172 2017-03-08 15:55:25Z knoop
329! Bugfix: moved parallel random generator initialization into its module
330!
331! 2118 2017-01-17 16:38:49Z raasch
332! OpenACC directives removed
333!
334! 2037 2016-10-26 11:15:40Z knoop
335! Anelastic approximation implemented
336!
337! 2031 2016-10-21 15:11:58Z knoop
338! renamed variable rho to rho_ocean
339!
340! 2011 2016-09-19 17:29:57Z kanani
341! Flag urban_surface is now defined in module control_parameters.
342!
343! 2007 2016-08-24 15:47:17Z kanani
344! Added support for urban surface model,
345! adjusted location_message in case of plant_canopy
346!
347! 2000 2016-08-20 18:09:15Z knoop
348! Forced header and separation lines into 80 columns
349!
350! 1992 2016-08-12 15:14:59Z suehring
351! Initializaton of scalarflux at model top
352! Bugfixes in initialization of surface and top salinity flux, top scalar and
353! humidity fluxes
354!
355! 1960 2016-07-12 16:34:24Z suehring
356! Separate humidity and passive scalar
357! Increase dimension for mean_inflow_profiles
358! Remove inadvertent write-statement
359! Bugfix, large-scale forcing is still not implemented for passive scalars
360!
361! 1957 2016-07-07 10:43:48Z suehring
362! flight module added
363!
364! 1920 2016-05-30 10:50:15Z suehring
365! Initialize us with very small number to avoid segmentation fault during
366! calculation of Obukhov length
367!
368! 1918 2016-05-27 14:35:57Z raasch
369! intermediate_timestep_count is set 0 instead 1 for first call of pres,
370! bugfix: initialization of local sum arrays are moved to the beginning of the
371!         routine because otherwise results from pres are overwritten
372!
373! 1914 2016-05-26 14:44:07Z witha
374! Added initialization of the wind turbine model
375!
376! 1878 2016-04-19 12:30:36Z hellstea
377! The zeroth element of weight_pres removed as unnecessary
378!
379! 1849 2016-04-08 11:33:18Z hoffmann
380! Adapted for modularization of microphysics.
381! precipitation_amount, precipitation_rate, prr moved to arrays_3d.
382! Initialization of nc_1d, nr_1d, pt_1d, qc_1d, qr_1d, q_1d moved to
383! bcm_init.
384!
385! 1845 2016-04-08 08:29:13Z raasch
386! nzb_2d replaced by nzb_u|v_inner
387!
388! 1833 2016-04-07 14:23:03Z raasch
389! initialization of spectra quantities moved to spectra_mod
390!
391! 1831 2016-04-07 13:15:51Z hoffmann
392! turbulence renamed collision_turbulence
393!
394! 1826 2016-04-07 12:01:39Z maronga
395! Renamed radiation calls.
396! Renamed canopy model calls.
397!
398! 1822 2016-04-07 07:49:42Z hoffmann
399! icloud_scheme replaced by microphysics_*
400!
401! 1817 2016-04-06 15:44:20Z maronga
402! Renamed lsm calls.
403!
404! 1815 2016-04-06 13:49:59Z raasch
405! zero-settings for velocities inside topography re-activated (was deactivated
406! in r1762)
407!
408! 1788 2016-03-10 11:01:04Z maronga
409! Added z0q.
410! Syntax layout improved.
411!
412! 1783 2016-03-06 18:36:17Z raasch
413! netcdf module name changed + related changes
414!
415! 1764 2016-02-28 12:45:19Z raasch
416! bugfix: increase size of volume_flow_area_l and volume_flow_initial_l by 1
417!
418! 1762 2016-02-25 12:31:13Z hellstea
419! Introduction of nested domain feature
420!
421! 1738 2015-12-18 13:56:05Z raasch
422! calculate mean surface level height for each statistic region
423!
424! 1734 2015-12-02 12:17:12Z raasch
425! no initial disturbances in case that the disturbance energy limit has been
426! set zero
427!
428! 1707 2015-11-02 15:24:52Z maronga
429! Bugfix: transfer of Richardson number from 1D model to Obukhov length caused
430! devision by zero in neutral stratification
431!
432! 1691 2015-10-26 16:17:44Z maronga
433! Call to init_surface_layer added. rif is replaced by ol and zeta.
434!
435! 1682 2015-10-07 23:56:08Z knoop
436! Code annotations made doxygen readable
437!
438! 1615 2015-07-08 18:49:19Z suehring
439! Enable turbulent inflow for passive_scalar and humidity
440!
441! 1585 2015-04-30 07:05:52Z maronga
442! Initialization of radiation code is now done after LSM initializtion
443!
444! 1575 2015-03-27 09:56:27Z raasch
445! adjustments for psolver-queries
446!
447! 1551 2015-03-03 14:18:16Z maronga
448! Allocation of land surface arrays is now done in the subroutine lsm_init_arrays,
449! which is part of land_surface_model.
450!
451! 1507 2014-12-10 12:14:18Z suehring
452! Bugfix: set horizontal velocity components to zero inside topography
453!
454! 1496 2014-12-02 17:25:50Z maronga
455! Added initialization of the land surface and radiation schemes
456!
457! 1484 2014-10-21 10:53:05Z kanani
458! Changes due to new module structure of the plant canopy model:
459! canopy-related initialization (e.g. lad and canopy_heat_flux) moved to new
460! subroutine init_plant_canopy within the module plant_canopy_model_mod,
461! call of subroutine init_plant_canopy added.
462!
463! 1431 2014-07-15 14:47:17Z suehring
464! var_d added, in order to normalize spectra.
465!
466! 1429 2014-07-15 12:53:45Z knoop
467! Ensemble run capability added to parallel random number generator
468!
469! 1411 2014-05-16 18:01:51Z suehring
470! Initial horizontal velocity profiles were not set to zero at the first vertical
471! grid level in case of non-cyclic lateral boundary conditions.
472!
473! 1406 2014-05-16 13:47:01Z raasch
474! bugfix: setting of initial velocities at k=1 to zero not in case of a
475! no-slip boundary condition for uv
476!
477! 1402 2014-05-09 14:25:13Z raasch
478! location messages modified
479!
480! 1400 2014-05-09 14:03:54Z knoop
481! Parallel random number generator added
482!
483! 1384 2014-05-02 14:31:06Z raasch
484! location messages added
485!
486! 1361 2014-04-16 15:17:48Z hoffmann
487! tend_* removed
488! Bugfix: w_subs is not allocated anymore if it is already allocated
489!
490! 1359 2014-04-11 17:15:14Z hoffmann
491! module lpm_init_mod added to use statements, because lpm_init has become a
492! module
493!
494! 1353 2014-04-08 15:21:23Z heinze
495! REAL constants provided with KIND-attribute
496!
497! 1340 2014-03-25 19:45:13Z kanani
498! REAL constants defined as wp-kind
499!
500! 1322 2014-03-20 16:38:49Z raasch
501! REAL constants defined as wp-kind
502! module interfaces removed
503!
504! 1320 2014-03-20 08:40:49Z raasch
505! ONLY-attribute added to USE-statements,
506! kind-parameters added to all INTEGER and REAL declaration statements,
507! kinds are defined in new module kinds,
508! revision history before 2012 removed,
509! comment fields (!:) to be used for variable explanations added to
510! all variable declaration statements
511!
512! 1316 2014-03-17 07:44:59Z heinze
513! Bugfix: allocation of w_subs
514!
515! 1299 2014-03-06 13:15:21Z heinze
516! Allocate w_subs due to extension of large scale subsidence in combination
517! with large scale forcing data (LSF_DATA)
518!
519! 1241 2013-10-30 11:36:58Z heinze
520! Overwrite initial profiles in case of nudging
521! Inititialize shf and qsws in case of large_scale_forcing
522!
523! 1221 2013-09-10 08:59:13Z raasch
524! +rflags_s_inner in copyin statement, use copyin for most arrays instead of
525! copy
526!
527! 1212 2013-08-15 08:46:27Z raasch
528! array tri is allocated and included in data copy statement
529!
530! 1195 2013-07-01 12:27:57Z heinze
531! Bugfix: move allocation of ref_state to parin.f90 and read_var_list.f90
532!
533! 1179 2013-06-14 05:57:58Z raasch
534! allocate and set ref_state to be used in buoyancy terms
535!
536! 1171 2013-05-30 11:27:45Z raasch
537! diss array is allocated with full size if accelerator boards are used
538!
539! 1159 2013-05-21 11:58:22Z fricke
540! -bc_lr_dirneu, bc_lr_neudir, bc_ns_dirneu, bc_ns_neudir
541!
542! 1153 2013-05-10 14:33:08Z raasch
543! diss array is allocated with dummy elements even if it is not needed
544! (required by PGI 13.4 / CUDA 5.0)
545!
546! 1115 2013-03-26 18:16:16Z hoffmann
547! unused variables removed
548!
549! 1113 2013-03-10 02:48:14Z raasch
550! openACC directive modified
551!
552! 1111 2013-03-08 23:54:10Z raasch
553! openACC directives added for pres
554! array diss allocated only if required
555!
556! 1092 2013-02-02 11:24:22Z raasch
557! unused variables removed
558!
559! 1065 2012-11-22 17:42:36Z hoffmann
560! allocation of diss (dissipation rate) in case of turbulence = .TRUE. added
561!
562! 1053 2012-11-13 17:11:03Z hoffmann
563! allocation and initialisation of necessary data arrays for the two-moment
564! cloud physics scheme the two new prognostic equations (nr, qr):
565! +dr, lambda_r, mu_r, sed_*, xr, *s, *sws, *swst, *, *_p, t*_m, *_1, *_2, *_3,
566! +tend_*, prr
567!
568! 1036 2012-10-22 13:43:42Z raasch
569! code put under GPL (PALM 3.9)
570!
571! 1032 2012-10-21 13:03:21Z letzel
572! save memory by not allocating pt_2 in case of neutral = .T.
573!
574! 1025 2012-10-07 16:04:41Z letzel
575! bugfix: swap indices of mask for ghost boundaries
576!
577! 1015 2012-09-27 09:23:24Z raasch
578! mask is set to zero for ghost boundaries
579!
580! 1010 2012-09-20 07:59:54Z raasch
581! cpp switch __nopointer added for pointer free version
582!
583! 1003 2012-09-14 14:35:53Z raasch
584! nxra,nyna, nzta replaced ny nxr, nyn, nzt
585!
586! 1001 2012-09-13 14:08:46Z raasch
587! all actions concerning leapfrog scheme removed
588!
589! 996 2012-09-07 10:41:47Z raasch
590! little reformatting
591!
592! 978 2012-08-09 08:28:32Z fricke
593! outflow damping layer removed
594! roughness length for scalar quantites z0h added
595! damping zone for the potential temperatur in case of non-cyclic lateral
596! boundaries added
597! initialization of ptdf_x, ptdf_y
598! initialization of c_u_m, c_u_m_l, c_v_m, c_v_m_l, c_w_m, c_w_m_l
599!
600! 849 2012-03-15 10:35:09Z raasch
601! init_particles renamed lpm_init
602!
603! 825 2012-02-19 03:03:44Z raasch
604! wang_collision_kernel renamed wang_kernel
605!
606! Revision 1.1  1998/03/09 16:22:22  raasch
607! Initial revision
608!
609!
610! Description:
611! ------------
612!> Allocation of arrays and initialization of the 3D model via
613!> a) pre-run the 1D model
614!> or
615!> b) pre-set constant linear profiles
616!> or
617!> c) read values of a previous run
618!------------------------------------------------------------------------------!
619 SUBROUTINE init_3d_model
620
621
622    USE advec_ws
623
624    USE arrays_3d
625
626    USE basic_constants_and_equations_mod,                                     &
627        ONLY:  c_p, g, l_v, pi, exner_function, exner_function_invers,         &
628               ideal_gas_law_rho, ideal_gas_law_rho_pt, barometric_formula
629
630    USE bulk_cloud_model_mod,                                                  &
631        ONLY:  bulk_cloud_model
632
633    USE chem_modules,                                                          &
634        ONLY:  max_pr_cs ! ToDo: this dependency needs to be removed cause it is ugly #new_dom
635
636    USE control_parameters
637
638    USE grid_variables,                                                        &
639        ONLY:  dx, dy, ddx2_mg, ddy2_mg
640
641    USE indices
642
643    USE kinds
644 
645    USE lsf_nudging_mod,                                                       &
646        ONLY:  ls_forcing_surf
647
648    USE model_1d_mod,                                                          &
649        ONLY:  init_1d_model, l1d, u1d, v1d
650
651    USE module_interface,                                                      &
652        ONLY:  module_interface_init_arrays,                                   &
653               module_interface_init,                                          &
654               module_interface_init_checks
655
656    USE multi_agent_system_mod,                                                &
657        ONLY:  agents_active, mas_init
658
659    USE netcdf_interface,                                                      &
660        ONLY:  dots_max
661
662    USE netcdf_data_input_mod,                                                 &
663        ONLY:  init_3d,                              &
664               netcdf_data_input_init_3d, netcdf_data_input_interpolate
665
666    USE nesting_offl_mod,                                                      &
667        ONLY:  nesting_offl_init
668
669    USE pegrid
670
671#if defined( __parallel )
672    USE pmc_interface,                                                         &
673        ONLY:  nested_run
674#endif
675
676    USE random_function_mod
677
678    USE random_generator_parallel,                                             &
679        ONLY:  init_parallel_random_generator
680
681    USE read_restart_data_mod,                                                 &
682        ONLY:  rrd_read_parts_of_global, rrd_local
683
684    USE statistics,                                                            &
685        ONLY:  hom, hom_sum, mean_surface_level_height, pr_palm, rmask,        &
686               statistic_regions, sums, sums_divnew_l, sums_divold_l, sums_l,  &
687               sums_l_l, sums_wsts_bc_l, ts_value,                             &
688               weight_pres, weight_substep
689
690    USE synthetic_turbulence_generator_mod,                                    &
691        ONLY:  stg_init, use_syn_turb_gen
692
693    USE surface_layer_fluxes_mod,                                              &
694        ONLY:  init_surface_layer_fluxes
695
696    USE surface_mod,                                                           &
697        ONLY :  init_surface_arrays, init_surfaces, surf_def_h, surf_lsm_h,    &
698                surf_usm_h, get_topography_top_index_ji
699
700#if defined( _OPENACC )
701    USE surface_mod,                                                           &
702        ONLY :  bc_h
703#endif
704
705    USE surface_data_output_mod,                                               &
706        ONLY:  surface_data_output_init
707
708    USE transpose_indices
709
710    IMPLICIT NONE
711
712    INTEGER(iwp) ::  i             !< grid index in x direction
713    INTEGER(iwp) ::  ind_array(1)  !< dummy used to determine start index for external pressure forcing
714    INTEGER(iwp) ::  j             !< grid index in y direction
715    INTEGER(iwp) ::  k             !< grid index in z direction
716    INTEGER(iwp) ::  k_surf        !< surface level index
717    INTEGER(iwp) ::  m             !< index of surface element in surface data type
718    INTEGER(iwp) ::  sr            !< index of statistic region
719
720    INTEGER(iwp), DIMENSION(:), ALLOCATABLE   ::  ngp_2dh_l  !< toal number of horizontal grid points in statistical region on subdomain
721
722    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  ngp_2dh_outer_l    !< number of horizontal non-wall bounded grid points on subdomain
723    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  ngp_2dh_s_inner_l  !< number of horizontal non-topography grid points on subdomain
724
725    REAL(wp), DIMENSION(:), ALLOCATABLE ::  init_l        !< dummy array used for averaging 3D data to obtain inital profiles
726    REAL(wp), DIMENSION(:), ALLOCATABLE ::  p_hydrostatic !< hydrostatic pressure
727
728    INTEGER(iwp) ::  l       !< loop variable
729    INTEGER(iwp) ::  nzt_l   !< index of top PE boundary for multigrid level
730    REAL(wp) ::  dx_l !< grid spacing along x on different multigrid level
731    REAL(wp) ::  dy_l !< grid spacing along y on different multigrid level
732
733    REAL(wp), DIMENSION(1:3) ::  volume_flow_area_l     !< area of lateral and top model domain surface on local subdomain
734    REAL(wp), DIMENSION(1:3) ::  volume_flow_initial_l  !< initial volume flow into model domain
735
736    REAL(wp), DIMENSION(:), ALLOCATABLE ::  mean_surface_level_height_l !< mean surface level height on subdomain
737    REAL(wp), DIMENSION(:), ALLOCATABLE ::  ngp_3d_inner_l    !< total number of non-topography grid points on subdomain
738    REAL(wp), DIMENSION(:), ALLOCATABLE ::  ngp_3d_inner_tmp  !< total number of non-topography grid points
739
740    INTEGER(iwp) ::  nz_u_shift   !< topography-top index on u-grid, used to vertically shift initial profiles
741    INTEGER(iwp) ::  nz_v_shift   !< topography-top index on v-grid, used to vertically shift initial profiles
742    INTEGER(iwp) ::  nz_w_shift   !< topography-top index on w-grid, used to vertically shift initial profiles
743    INTEGER(iwp) ::  nz_s_shift   !< topography-top index on scalar-grid, used to vertically shift initial profiles
744    INTEGER(iwp) ::  nz_u_shift_l !< topography-top index on u-grid, used to vertically shift initial profiles
745    INTEGER(iwp) ::  nz_v_shift_l !< topography-top index on v-grid, used to vertically shift initial profiles
746    INTEGER(iwp) ::  nz_w_shift_l !< topography-top index on w-grid, used to vertically shift initial profiles
747    INTEGER(iwp) ::  nz_s_shift_l !< topography-top index on scalar-grid, used to vertically shift initial profiles
748
749
750    CALL location_message( 'model initialization', 'start' )
751
752    IF ( debug_output )  CALL debug_message( 'allocating arrays', 'start' )
753!
754!-- Allocate arrays
755    ALLOCATE( mean_surface_level_height(0:statistic_regions),                  &
756              mean_surface_level_height_l(0:statistic_regions),                &
757              ngp_2dh(0:statistic_regions), ngp_2dh_l(0:statistic_regions),    &
758              ngp_3d(0:statistic_regions),                                     &
759              ngp_3d_inner(0:statistic_regions),                               &
760              ngp_3d_inner_l(0:statistic_regions),                             &
761              ngp_3d_inner_tmp(0:statistic_regions),                           &
762              sums_divnew_l(0:statistic_regions),                              &
763              sums_divold_l(0:statistic_regions) )
764    ALLOCATE( dp_smooth_factor(nzb:nzt), rdf(nzb+1:nzt), rdf_sc(nzb+1:nzt) )
765    ALLOCATE( ngp_2dh_outer(nzb:nzt+1,0:statistic_regions),                    &
766              ngp_2dh_outer_l(nzb:nzt+1,0:statistic_regions),                  &
767              ngp_2dh_s_inner(nzb:nzt+1,0:statistic_regions),                  &
768              ngp_2dh_s_inner_l(nzb:nzt+1,0:statistic_regions),                &
769              rmask(nysg:nyng,nxlg:nxrg,0:statistic_regions),                  &
770              sums(nzb:nzt+1,pr_palm+max_pr_user+max_pr_cs+max_pr_salsa),      &
771              sums_l(nzb:nzt+1,pr_palm+max_pr_user+max_pr_cs+max_pr_salsa,0:threads_per_task-1),   &
772              sums_l_l(nzb:nzt+1,0:statistic_regions,0:threads_per_task-1),    &
773              sums_wsts_bc_l(nzb:nzt+1,0:statistic_regions) )
774    ALLOCATE( ts_value(dots_max,0:statistic_regions) )
775    ALLOCATE( ptdf_x(nxlg:nxrg), ptdf_y(nysg:nyng) )
776
777    ALLOCATE( d(nzb+1:nzt,nys:nyn,nxl:nxr),                                    &
778              p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                &
779              tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
780
781    ALLOCATE( pt_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
782              pt_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
783              u_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
784              u_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
785              u_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
786              v_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
787              v_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
788              v_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
789              w_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
790              w_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
791              w_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
792    IF (  .NOT.  neutral )  THEN
793       ALLOCATE( pt_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
794    ENDIF
795!
796!-- Pre-set masks for regional statistics. Default is the total model domain.
797!-- Ghost points are excluded because counting values at the ghost boundaries
798!-- would bias the statistics
799    rmask = 1.0_wp
800    rmask(:,nxlg:nxl-1,:) = 0.0_wp;  rmask(:,nxr+1:nxrg,:) = 0.0_wp
801    rmask(nysg:nys-1,:,:) = 0.0_wp;  rmask(nyn+1:nyng,:,:) = 0.0_wp
802!
803!-- Following array is required for perturbation pressure within the iterative
804!-- pressure solvers. For the multistep schemes (Runge-Kutta), array p holds
805!-- the weighted average of the substeps and cannot be used in the Poisson
806!-- solver.
807    IF ( psolver == 'sor' )  THEN
808       ALLOCATE( p_loc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
809    ELSEIF ( psolver(1:9) == 'multigrid' )  THEN
810!
811!--    For performance reasons, multigrid is using one ghost layer only
812       ALLOCATE( p_loc(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
813    ENDIF
814
815!
816!-- Array for storing constant coeffficients of the tridiagonal solver
817    IF ( psolver == 'poisfft' )  THEN
818       ALLOCATE( tri(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1,2) )
819       ALLOCATE( tric(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1) )
820    ENDIF
821
822    IF ( humidity )  THEN
823!
824!--    3D-humidity
825       ALLOCATE( q_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
826                 q_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
827                 q_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
828                 vpt_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 
829    ENDIF   
830   
831    IF ( passive_scalar )  THEN
832
833!
834!--    3D scalar arrays
835       ALLOCATE( s_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
836                 s_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
837                 s_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
838
839    ENDIF
840
841!
842!-- Allocate and set 1d-profiles for Stokes drift velocity. It may be set to
843!-- non-zero values later in ocean_init
844    ALLOCATE( u_stokes_zu(nzb:nzt+1), u_stokes_zw(nzb:nzt+1),                  &
845              v_stokes_zu(nzb:nzt+1), v_stokes_zw(nzb:nzt+1) )
846    u_stokes_zu(:) = 0.0_wp
847    u_stokes_zw(:) = 0.0_wp
848    v_stokes_zu(:) = 0.0_wp
849    v_stokes_zw(:) = 0.0_wp
850
851!
852!-- Allocation of anelastic and Boussinesq approximation specific arrays
853    ALLOCATE( p_hydrostatic(nzb:nzt+1) )
854    ALLOCATE( rho_air(nzb:nzt+1) )
855    ALLOCATE( rho_air_zw(nzb:nzt+1) )
856    ALLOCATE( drho_air(nzb:nzt+1) )
857    ALLOCATE( drho_air_zw(nzb:nzt+1) )
858!
859!-- Density profile calculation for anelastic and Boussinesq approximation
860!-- In case of a Boussinesq approximation, a constant density is calculated
861!-- mainly for output purposes. This density do not need to be considered
862!-- in the model's system of equations.
863    IF ( TRIM( approximation ) == 'anelastic' )  THEN
864       DO  k = nzb, nzt+1
865          p_hydrostatic(k) = barometric_formula(zu(k), pt_surface *            & 
866                             exner_function(surface_pressure * 100.0_wp),      &
867                             surface_pressure * 100.0_wp)
868         
869          rho_air(k) = ideal_gas_law_rho_pt(p_hydrostatic(k), pt_init(k))
870       ENDDO
871       
872       DO  k = nzb, nzt
873          rho_air_zw(k) = 0.5_wp * ( rho_air(k) + rho_air(k+1) )
874       ENDDO
875       
876       rho_air_zw(nzt+1)  = rho_air_zw(nzt)                                    &
877                            + 2.0_wp * ( rho_air(nzt+1) - rho_air_zw(nzt)  )
878                           
879    ELSE
880       DO  k = nzb, nzt+1
881          p_hydrostatic(k) = barometric_formula(zu(nzb), pt_surface *          & 
882                             exner_function(surface_pressure * 100.0_wp),      &
883                             surface_pressure * 100.0_wp)
884
885          rho_air(k) = ideal_gas_law_rho_pt(p_hydrostatic(k), pt_init(nzb))
886       ENDDO
887       
888       DO  k = nzb, nzt
889          rho_air_zw(k) = 0.5_wp * ( rho_air(k) + rho_air(k+1) )
890       ENDDO
891       
892       rho_air_zw(nzt+1)  = rho_air_zw(nzt)                                    &
893                            + 2.0_wp * ( rho_air(nzt+1) - rho_air_zw(nzt)  )
894                           
895    ENDIF
896!
897!-- compute the inverse density array in order to avoid expencive divisions
898    drho_air    = 1.0_wp / rho_air
899    drho_air_zw = 1.0_wp / rho_air_zw
900
901!
902!-- Allocation of flux conversion arrays
903    ALLOCATE( heatflux_input_conversion(nzb:nzt+1) )
904    ALLOCATE( waterflux_input_conversion(nzb:nzt+1) )
905    ALLOCATE( momentumflux_input_conversion(nzb:nzt+1) )
906    ALLOCATE( heatflux_output_conversion(nzb:nzt+1) )
907    ALLOCATE( waterflux_output_conversion(nzb:nzt+1) )
908    ALLOCATE( momentumflux_output_conversion(nzb:nzt+1) )
909
910!
911!-- calculate flux conversion factors according to approximation and in-/output mode
912    DO  k = nzb, nzt+1
913
914        IF ( TRIM( flux_input_mode ) == 'kinematic' )  THEN
915            heatflux_input_conversion(k)      = rho_air_zw(k)
916            waterflux_input_conversion(k)     = rho_air_zw(k)
917            momentumflux_input_conversion(k)  = rho_air_zw(k)
918        ELSEIF ( TRIM( flux_input_mode ) == 'dynamic' ) THEN
919            heatflux_input_conversion(k)      = 1.0_wp / c_p
920            waterflux_input_conversion(k)     = 1.0_wp / l_v
921            momentumflux_input_conversion(k)  = 1.0_wp
922        ENDIF
923
924        IF ( TRIM( flux_output_mode ) == 'kinematic' )  THEN
925            heatflux_output_conversion(k)     = drho_air_zw(k)
926            waterflux_output_conversion(k)    = drho_air_zw(k)
927            momentumflux_output_conversion(k) = drho_air_zw(k)
928        ELSEIF ( TRIM( flux_output_mode ) == 'dynamic' ) THEN
929            heatflux_output_conversion(k)     = c_p
930            waterflux_output_conversion(k)    = l_v
931            momentumflux_output_conversion(k) = 1.0_wp
932        ENDIF
933
934        IF ( .NOT. humidity ) THEN
935            waterflux_input_conversion(k)  = 1.0_wp
936            waterflux_output_conversion(k) = 1.0_wp
937        ENDIF
938
939    ENDDO
940
941!
942!-- In case of multigrid method, compute grid lengths and grid factors for the
943!-- grid levels with respective density on each grid
944    IF ( psolver(1:9) == 'multigrid' )  THEN
945
946       ALLOCATE( ddx2_mg(maximum_grid_level) )
947       ALLOCATE( ddy2_mg(maximum_grid_level) )
948       ALLOCATE( dzu_mg(nzb+1:nzt+1,maximum_grid_level) )
949       ALLOCATE( dzw_mg(nzb+1:nzt+1,maximum_grid_level) )
950       ALLOCATE( f1_mg(nzb+1:nzt,maximum_grid_level) )
951       ALLOCATE( f2_mg(nzb+1:nzt,maximum_grid_level) )
952       ALLOCATE( f3_mg(nzb+1:nzt,maximum_grid_level) )
953       ALLOCATE( rho_air_mg(nzb:nzt+1,maximum_grid_level) )
954       ALLOCATE( rho_air_zw_mg(nzb:nzt+1,maximum_grid_level) )
955
956       dzu_mg(:,maximum_grid_level) = dzu
957       rho_air_mg(:,maximum_grid_level) = rho_air
958!       
959!--    Next line to ensure an equally spaced grid.
960       dzu_mg(1,maximum_grid_level) = dzu(2)
961       rho_air_mg(nzb,maximum_grid_level) = rho_air(nzb) +                     &
962                                             (rho_air(nzb) - rho_air(nzb+1))
963
964       dzw_mg(:,maximum_grid_level) = dzw
965       rho_air_zw_mg(:,maximum_grid_level) = rho_air_zw
966       nzt_l = nzt
967       DO  l = maximum_grid_level-1, 1, -1
968           dzu_mg(nzb+1,l) = 2.0_wp * dzu_mg(nzb+1,l+1)
969           dzw_mg(nzb+1,l) = 2.0_wp * dzw_mg(nzb+1,l+1)
970           rho_air_mg(nzb,l)    = rho_air_mg(nzb,l+1) + (rho_air_mg(nzb,l+1) - rho_air_mg(nzb+1,l+1))
971           rho_air_zw_mg(nzb,l) = rho_air_zw_mg(nzb,l+1) + (rho_air_zw_mg(nzb,l+1) - rho_air_zw_mg(nzb+1,l+1))
972           rho_air_mg(nzb+1,l)    = rho_air_mg(nzb+1,l+1)
973           rho_air_zw_mg(nzb+1,l) = rho_air_zw_mg(nzb+1,l+1)
974           nzt_l = nzt_l / 2
975           DO  k = 2, nzt_l+1
976              dzu_mg(k,l) = dzu_mg(2*k-2,l+1) + dzu_mg(2*k-1,l+1)
977              dzw_mg(k,l) = dzw_mg(2*k-2,l+1) + dzw_mg(2*k-1,l+1)
978              rho_air_mg(k,l)    = rho_air_mg(2*k-1,l+1)
979              rho_air_zw_mg(k,l) = rho_air_zw_mg(2*k-1,l+1)
980           ENDDO
981       ENDDO
982
983       nzt_l = nzt
984       dx_l  = dx
985       dy_l  = dy
986       DO  l = maximum_grid_level, 1, -1
987          ddx2_mg(l) = 1.0_wp / dx_l**2
988          ddy2_mg(l) = 1.0_wp / dy_l**2
989          DO  k = nzb+1, nzt_l
990             f2_mg(k,l) = rho_air_zw_mg(k,l) / ( dzu_mg(k+1,l) * dzw_mg(k,l) )
991             f3_mg(k,l) = rho_air_zw_mg(k-1,l) / ( dzu_mg(k,l)   * dzw_mg(k,l) )
992             f1_mg(k,l) = 2.0_wp * ( ddx2_mg(l) + ddy2_mg(l) ) &
993                          * rho_air_mg(k,l) + f2_mg(k,l) + f3_mg(k,l)
994          ENDDO
995          nzt_l = nzt_l / 2
996          dx_l  = dx_l * 2.0_wp
997          dy_l  = dy_l * 2.0_wp
998       ENDDO
999
1000    ENDIF
1001
1002!
1003!-- 1D-array for large scale subsidence velocity
1004    IF ( .NOT. ALLOCATED( w_subs ) )  THEN
1005       ALLOCATE ( w_subs(nzb:nzt+1) )
1006       w_subs = 0.0_wp
1007    ENDIF
1008
1009!
1010!-- Arrays to store velocity data from t-dt and the phase speeds which
1011!-- are needed for radiation boundary conditions
1012    IF ( bc_radiation_l )  THEN
1013       ALLOCATE( u_m_l(nzb:nzt+1,nysg:nyng,1:2),                               &
1014                 v_m_l(nzb:nzt+1,nysg:nyng,0:1),                               &
1015                 w_m_l(nzb:nzt+1,nysg:nyng,0:1) )
1016    ENDIF
1017    IF ( bc_radiation_r )  THEN
1018       ALLOCATE( u_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx),                           &
1019                 v_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx),                           &
1020                 w_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx) )
1021    ENDIF
1022    IF ( bc_radiation_l  .OR.  bc_radiation_r )  THEN
1023       ALLOCATE( c_u(nzb:nzt+1,nysg:nyng), c_v(nzb:nzt+1,nysg:nyng),           &
1024                 c_w(nzb:nzt+1,nysg:nyng) )
1025    ENDIF
1026    IF ( bc_radiation_s )  THEN
1027       ALLOCATE( u_m_s(nzb:nzt+1,0:1,nxlg:nxrg),                               &
1028                 v_m_s(nzb:nzt+1,1:2,nxlg:nxrg),                               &
1029                 w_m_s(nzb:nzt+1,0:1,nxlg:nxrg) )
1030    ENDIF
1031    IF ( bc_radiation_n )  THEN
1032       ALLOCATE( u_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg),                           &
1033                 v_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg),                           &
1034                 w_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg) )
1035    ENDIF
1036    IF ( bc_radiation_s  .OR.  bc_radiation_n )  THEN
1037       ALLOCATE( c_u(nzb:nzt+1,nxlg:nxrg), c_v(nzb:nzt+1,nxlg:nxrg),           &
1038                 c_w(nzb:nzt+1,nxlg:nxrg) )
1039    ENDIF
1040    IF ( bc_radiation_l  .OR.  bc_radiation_r  .OR.  bc_radiation_s  .OR.      &
1041         bc_radiation_n )  THEN
1042       ALLOCATE( c_u_m_l(nzb:nzt+1), c_v_m_l(nzb:nzt+1), c_w_m_l(nzb:nzt+1) )                   
1043       ALLOCATE( c_u_m(nzb:nzt+1), c_v_m(nzb:nzt+1), c_w_m(nzb:nzt+1) )
1044    ENDIF
1045
1046!
1047!-- Initial assignment of the pointers
1048    IF ( .NOT. neutral )  THEN
1049       pt => pt_1;  pt_p => pt_2;  tpt_m => pt_3
1050    ELSE
1051       pt => pt_1;  pt_p => pt_1;  tpt_m => pt_3
1052    ENDIF
1053    u  => u_1;   u_p  => u_2;   tu_m  => u_3
1054    v  => v_1;   v_p  => v_2;   tv_m  => v_3
1055    w  => w_1;   w_p  => w_2;   tw_m  => w_3
1056
1057    IF ( humidity )  THEN
1058       q => q_1;  q_p => q_2;  tq_m => q_3
1059       vpt  => vpt_1
1060    ENDIF
1061   
1062    IF ( passive_scalar )  THEN
1063       s => s_1;  s_p => s_2;  ts_m => s_3
1064    ENDIF   
1065
1066!
1067!-- Initialize surface arrays
1068    CALL init_surface_arrays
1069!
1070!-- Allocate arrays for other modules
1071    CALL module_interface_init_arrays
1072
1073
1074!
1075!-- Allocate arrays containing the RK coefficient for calculation of
1076!-- perturbation pressure and turbulent fluxes. At this point values are
1077!-- set for pressure calculation during initialization (where no timestep
1078!-- is done). Further below the values needed within the timestep scheme
1079!-- will be set.
1080    ALLOCATE( weight_substep(1:intermediate_timestep_count_max),               &
1081              weight_pres(1:intermediate_timestep_count_max) )
1082    weight_substep = 1.0_wp
1083    weight_pres    = 1.0_wp
1084    intermediate_timestep_count = 0  ! needed when simulated_time = 0.0
1085       
1086    IF ( debug_output )  CALL debug_message( 'allocating arrays', 'end' )
1087
1088!
1089!-- Initialize time series
1090    ts_value = 0.0_wp
1091
1092!
1093!-- Initialize local summation arrays for routine flow_statistics.
1094!-- This is necessary because they may not yet have been initialized when they
1095!-- are called from flow_statistics (or - depending on the chosen model run -
1096!-- are never initialized)
1097    sums_divnew_l      = 0.0_wp
1098    sums_divold_l      = 0.0_wp
1099    sums_l_l           = 0.0_wp
1100    sums_wsts_bc_l     = 0.0_wp
1101   
1102!
1103!-- Initialize model variables
1104    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
1105         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
1106!
1107!--    Initialization with provided input data derived from larger-scale model
1108       IF ( INDEX( initializing_actions, 'inifor' ) /= 0 )  THEN
1109          IF ( debug_output )  CALL debug_message( 'initializing with INIFOR', 'start' )
1110!
1111!--       Read initial 1D profiles or 3D data from NetCDF file, depending
1112!--       on the provided level-of-detail.
1113!--       At the moment, only u, v, w, pt and q are provided.
1114          CALL netcdf_data_input_init_3d
1115!
1116!--       Please note, Inifor provides data from nzb+1 to nzt.
1117!--       Bottom and top boundary conditions for Inifor profiles are already
1118!--       set (just after reading), so that this is not necessary here.
1119!--       Depending on the provided level-of-detail, initial Inifor data is
1120!--       either stored on data type (lod=1), or directly on 3D arrays (lod=2).
1121!--       In order to obtain also initial profiles in case of lod=2 (which
1122!--       is required for e.g. damping), average over 3D data.
1123          IF( init_3d%lod_u == 1 )  THEN
1124             u_init = init_3d%u_init
1125          ELSEIF( init_3d%lod_u == 2 )  THEN
1126             ALLOCATE( init_l(nzb:nzt+1) ) 
1127             DO  k = nzb, nzt+1
1128                init_l(k) = SUM( u(k,nys:nyn,nxl:nxr) )
1129             ENDDO
1130             init_l = init_l / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
1131
1132#if defined( __parallel )
1133             CALL MPI_ALLREDUCE( init_l, u_init, nzt+1-nzb+1,                  &
1134                                 MPI_REAL, MPI_SUM, comm2d, ierr )
1135#else
1136             u_init = init_l
1137#endif
1138             DEALLOCATE( init_l )
1139
1140          ENDIF
1141           
1142          IF( init_3d%lod_v == 1 )  THEN 
1143             v_init = init_3d%v_init
1144          ELSEIF( init_3d%lod_v == 2 )  THEN
1145             ALLOCATE( init_l(nzb:nzt+1) ) 
1146             DO  k = nzb, nzt+1
1147                init_l(k) = SUM( v(k,nys:nyn,nxl:nxr) )
1148             ENDDO
1149             init_l = init_l / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
1150
1151#if defined( __parallel )
1152             CALL MPI_ALLREDUCE( init_l, v_init, nzt+1-nzb+1,                  &
1153                                 MPI_REAL, MPI_SUM, comm2d, ierr )
1154#else
1155             v_init = init_l
1156#endif
1157             DEALLOCATE( init_l )
1158          ENDIF
1159          IF( .NOT. neutral )  THEN
1160             IF( init_3d%lod_pt == 1 )  THEN
1161                pt_init = init_3d%pt_init
1162             ELSEIF( init_3d%lod_pt == 2 )  THEN
1163                ALLOCATE( init_l(nzb:nzt+1) ) 
1164                DO  k = nzb, nzt+1
1165                   init_l(k) = SUM( pt(k,nys:nyn,nxl:nxr) )
1166                ENDDO
1167                init_l = init_l / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
1168
1169#if defined( __parallel )
1170                CALL MPI_ALLREDUCE( init_l, pt_init, nzt+1-nzb+1,               &
1171                                    MPI_REAL, MPI_SUM, comm2d, ierr )
1172#else
1173                pt_init = init_l
1174#endif
1175                DEALLOCATE( init_l )
1176             ENDIF
1177          ENDIF
1178
1179
1180          IF( humidity )  THEN
1181             IF( init_3d%lod_q == 1 )  THEN
1182                q_init = init_3d%q_init
1183             ELSEIF( init_3d%lod_q == 2 )  THEN
1184                ALLOCATE( init_l(nzb:nzt+1) ) 
1185                DO  k = nzb, nzt+1
1186                   init_l(k) = SUM( q(k,nys:nyn,nxl:nxr) )
1187                ENDDO
1188                init_l = init_l / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
1189
1190#if defined( __parallel )
1191                CALL MPI_ALLREDUCE( init_l, q_init, nzt+1-nzb+1,               &
1192                                    MPI_REAL, MPI_SUM, comm2d, ierr )
1193#else
1194                q_init = init_l
1195#endif
1196                DEALLOCATE( init_l )
1197             ENDIF
1198          ENDIF
1199
1200!
1201!--       Write initial profiles onto 3D arrays.
1202!--       Work-around, 3D initialization of u,v,w creates artificial
1203!--       structures wich correlate with the processor grid. The reason
1204!--       for this is still unknown. To work-around this, 3D initialization
1205!--       will be effectively reduce to a 1D initialization where no such
1206!--       artificial structures appear.
1207          DO  i = nxlg, nxrg
1208             DO  j = nysg, nyng
1209                IF( init_3d%lod_u == 1  .OR.  init_3d%lod_u == 2 )             &
1210                   u(:,j,i) = u_init(:)
1211                IF( init_3d%lod_v == 1  .OR.  init_3d%lod_u == 2 )             &
1212                   v(:,j,i) = v_init(:)
1213                IF( .NOT. neutral  .AND.                                       &
1214                    ( init_3d%lod_pt == 1  .OR.  init_3d%lod_pt == 2 ) )       &
1215                   pt(:,j,i) = pt_init(:)
1216                IF( humidity  .AND.                                            &
1217                    ( init_3d%lod_q == 1  .OR.  init_3d%lod_q == 2 ) )         &
1218                   q(:,j,i) = q_init(:)
1219             ENDDO
1220          ENDDO
1221!
1222!--       Set geostrophic wind components. 
1223          IF ( init_3d%from_file_ug )  THEN
1224             ug(:) = init_3d%ug_init(:)
1225          ENDIF
1226          IF ( init_3d%from_file_vg )  THEN
1227             vg(:) = init_3d%vg_init(:)
1228          ENDIF
1229!
1230!--       Set bottom and top boundary condition for geostrophic wind
1231          ug(nzt+1) = ug(nzt)
1232          vg(nzt+1) = vg(nzt)
1233          ug(nzb)   = ug(nzb+1)
1234          vg(nzb)   = vg(nzb+1)
1235!
1236!--       Set inital w to 0
1237          w = 0.0_wp
1238
1239          IF ( passive_scalar )  THEN
1240             DO  i = nxlg, nxrg
1241                DO  j = nysg, nyng
1242                   s(:,j,i) = s_init
1243                ENDDO
1244             ENDDO
1245          ENDIF
1246
1247!
1248!--       Set velocity components at non-atmospheric / oceanic grid points to
1249!--       zero.
1250          u = MERGE( u, 0.0_wp, BTEST( wall_flags_0, 1 ) )
1251          v = MERGE( v, 0.0_wp, BTEST( wall_flags_0, 2 ) )
1252          w = MERGE( w, 0.0_wp, BTEST( wall_flags_0, 3 ) )
1253!
1254!--       Initialize surface variables, e.g. friction velocity, momentum
1255!--       fluxes, etc.
1256          CALL init_surfaces
1257
1258          IF ( debug_output )  CALL debug_message( 'initializing with INIFOR', 'end' )
1259!
1260!--    Initialization via computed 1D-model profiles
1261       ELSEIF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
1262
1263          IF ( debug_output )  CALL debug_message( 'initializing with 1D model profiles', 'start' )
1264!
1265!--       Use solutions of the 1D model as initial profiles,
1266!--       start 1D model
1267          CALL init_1d_model
1268!
1269!--       Transfer initial profiles to the arrays of the 3D model
1270          DO  i = nxlg, nxrg
1271             DO  j = nysg, nyng
1272                pt(:,j,i) = pt_init
1273                u(:,j,i)  = u1d
1274                v(:,j,i)  = v1d
1275             ENDDO
1276          ENDDO
1277
1278          IF ( humidity )  THEN
1279             DO  i = nxlg, nxrg
1280                DO  j = nysg, nyng
1281                   q(:,j,i) = q_init
1282                ENDDO
1283             ENDDO
1284          ENDIF
1285
1286          IF ( passive_scalar )  THEN
1287             DO  i = nxlg, nxrg
1288                DO  j = nysg, nyng
1289                   s(:,j,i) = s_init
1290                ENDDO
1291             ENDDO   
1292          ENDIF
1293!
1294!--          Store initial profiles for output purposes etc.
1295          IF ( .NOT. constant_diffusion )  THEN
1296             hom(:,1,25,:) = SPREAD( l1d, 2, statistic_regions+1 )
1297          ENDIF
1298!
1299!--       Set velocities back to zero
1300          u = MERGE( u, 0.0_wp, BTEST( wall_flags_0, 1 ) )
1301          v = MERGE( v, 0.0_wp, BTEST( wall_flags_0, 2 ) )         
1302!
1303!--       WARNING: The extra boundary conditions set after running the
1304!--       -------  1D model impose an error on the divergence one layer
1305!--                below the topography; need to correct later
1306!--       ATTENTION: Provisional correction for Piacsek & Williams
1307!--       ---------  advection scheme: keep u and v zero one layer below
1308!--                  the topography.
1309          IF ( ibc_uv_b == 1 )  THEN
1310!
1311!--          Neumann condition
1312             DO  i = nxl-1, nxr+1
1313                DO  j = nys-1, nyn+1
1314                   u(nzb,j,i) = u(nzb+1,j,i)
1315                   v(nzb,j,i) = v(nzb+1,j,i)
1316                ENDDO
1317             ENDDO
1318
1319          ENDIF
1320!
1321!--       Initialize surface variables, e.g. friction velocity, momentum
1322!--       fluxes, etc.
1323          CALL init_surfaces
1324
1325          IF ( debug_output )  CALL debug_message( 'initializing with 1D model profiles', 'end' )
1326
1327       ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 )    &
1328       THEN
1329
1330          IF ( debug_output )  CALL debug_message( 'initializing with constant profiles', 'start' )
1331
1332!
1333!--       Use constructed initial profiles (velocity constant with height,
1334!--       temperature profile with constant gradient)
1335          DO  i = nxlg, nxrg
1336             DO  j = nysg, nyng
1337                pt(:,j,i) = pt_init
1338                u(:,j,i)  = u_init
1339                v(:,j,i)  = v_init
1340             ENDDO
1341          ENDDO
1342!
1343!--       Mask topography
1344          u = MERGE( u, 0.0_wp, BTEST( wall_flags_0, 1 ) )
1345          v = MERGE( v, 0.0_wp, BTEST( wall_flags_0, 2 ) )
1346!
1347!--       Set initial horizontal velocities at the lowest computational grid
1348!--       levels to zero in order to avoid too small time steps caused by the
1349!--       diffusion limit in the initial phase of a run (at k=1, dz/2 occurs
1350!--       in the limiting formula!).
1351!--       Please note, in case land- or urban-surface model is used and a
1352!--       spinup is applied, masking the lowest computational level is not
1353!--       possible as MOST as well as energy-balance parametrizations will not
1354!--       work with zero wind velocity.
1355          IF ( ibc_uv_b /= 1  .AND.  .NOT.  spinup )  THEN
1356             DO  i = nxlg, nxrg
1357                DO  j = nysg, nyng
1358                   DO  k = nzb, nzt
1359                      u(k,j,i) = MERGE( u(k,j,i), 0.0_wp,                      &
1360                                        BTEST( wall_flags_0(k,j,i), 20 ) )
1361                      v(k,j,i) = MERGE( v(k,j,i), 0.0_wp,                      &
1362                                        BTEST( wall_flags_0(k,j,i), 21 ) )
1363                   ENDDO
1364                ENDDO
1365             ENDDO
1366          ENDIF
1367
1368          IF ( humidity )  THEN
1369             DO  i = nxlg, nxrg
1370                DO  j = nysg, nyng
1371                   q(:,j,i) = q_init
1372                ENDDO
1373             ENDDO
1374          ENDIF
1375         
1376          IF ( passive_scalar )  THEN
1377             DO  i = nxlg, nxrg
1378                DO  j = nysg, nyng
1379                   s(:,j,i) = s_init
1380                ENDDO
1381             ENDDO
1382          ENDIF
1383
1384!
1385!--       Compute initial temperature field and other constants used in case
1386!--       of a sloping surface
1387          IF ( sloping_surface )  CALL init_slope
1388!
1389!--       Initialize surface variables, e.g. friction velocity, momentum
1390!--       fluxes, etc.
1391          CALL init_surfaces
1392         
1393          IF ( debug_output )  CALL debug_message( 'initializing with constant profiles', 'end' )
1394
1395       ELSEIF ( INDEX(initializing_actions, 'by_user') /= 0 )                  &
1396       THEN
1397
1398          IF ( debug_output )  CALL debug_message( 'initializing by user', 'start' )
1399!
1400!--       Pre-initialize surface variables, i.e. setting start- and end-indices
1401!--       at each (j,i)-location. Please note, this does not supersede
1402!--       user-defined initialization of surface quantities.
1403          CALL init_surfaces
1404!
1405!--       Initialization will completely be done by the user
1406          CALL user_init_3d_model
1407
1408          IF ( debug_output )  CALL debug_message( 'initializing by user', 'end' )
1409
1410       ENDIF
1411
1412       IF ( debug_output )  CALL debug_message( 'initializing statistics, boundary conditions, etc.', 'start' )
1413
1414!
1415!--    Bottom boundary
1416       IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2  )  THEN
1417          u(nzb,:,:) = 0.0_wp
1418          v(nzb,:,:) = 0.0_wp
1419       ENDIF
1420
1421!
1422!--    Apply channel flow boundary condition
1423       IF ( TRIM( bc_uv_t ) == 'dirichlet_0' )  THEN
1424          u(nzt+1,:,:) = 0.0_wp
1425          v(nzt+1,:,:) = 0.0_wp
1426       ENDIF
1427
1428!
1429!--    Calculate virtual potential temperature
1430       IF ( humidity )  vpt = pt * ( 1.0_wp + 0.61_wp * q )
1431
1432!
1433!--    Store initial profiles for output purposes etc.. Please note, in case of
1434!--    initialization of u, v, w, pt, and q via output data derived from larger
1435!--    scale models, data will not be horizontally homogeneous. Actually, a mean
1436!--    profile should be calculated before.   
1437       hom(:,1,5,:) = SPREAD( u(:,nys,nxl), 2, statistic_regions+1 )
1438       hom(:,1,6,:) = SPREAD( v(:,nys,nxl), 2, statistic_regions+1 )
1439       IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2)  THEN
1440          hom(nzb,1,5,:) = 0.0_wp
1441          hom(nzb,1,6,:) = 0.0_wp
1442       ENDIF
1443       hom(:,1,7,:)  = SPREAD( pt(:,nys,nxl), 2, statistic_regions+1 )
1444
1445       IF ( humidity )  THEN
1446!
1447!--       Store initial profile of total water content, virtual potential
1448!--       temperature
1449          hom(:,1,26,:) = SPREAD(   q(:,nys,nxl), 2, statistic_regions+1 )
1450          hom(:,1,29,:) = SPREAD( vpt(:,nys,nxl), 2, statistic_regions+1 )
1451!
1452!--       Store initial profile of mixing ratio and potential
1453!--       temperature
1454          IF ( bulk_cloud_model  .OR.  cloud_droplets ) THEN
1455             hom(:,1,27,:) = SPREAD(  q(:,nys,nxl), 2, statistic_regions+1 )
1456             hom(:,1,28,:) = SPREAD( pt(:,nys,nxl), 2, statistic_regions+1 )
1457          ENDIF
1458       ENDIF
1459
1460!
1461!--    Store initial scalar profile
1462       IF ( passive_scalar )  THEN
1463          hom(:,1,121,:) = SPREAD(  s(:,nys,nxl), 2, statistic_regions+1 )
1464       ENDIF
1465
1466!
1467!--    Initialize the random number generators (from numerical recipes)
1468       CALL random_function_ini
1469       
1470       IF ( random_generator == 'random-parallel' )  THEN
1471          CALL init_parallel_random_generator( nx, nys, nyn, nxl, nxr )
1472       ENDIF
1473!
1474!--    Set the reference state to be used in the buoyancy terms (for ocean runs
1475!--    the reference state will be set (overwritten) in init_ocean)
1476       IF ( use_single_reference_value )  THEN
1477          IF (  .NOT.  humidity )  THEN
1478             ref_state(:) = pt_reference
1479          ELSE
1480             ref_state(:) = vpt_reference
1481          ENDIF
1482       ELSE
1483          IF (  .NOT.  humidity )  THEN
1484             ref_state(:) = pt_init(:)
1485          ELSE
1486             ref_state(:) = vpt(:,nys,nxl)
1487          ENDIF
1488       ENDIF
1489
1490!
1491!--    For the moment, vertical velocity is zero
1492       w = 0.0_wp
1493
1494!
1495!--    Initialize array sums (must be defined in first call of pres)
1496       sums = 0.0_wp
1497
1498!
1499!--    In case of iterative solvers, p must get an initial value
1500       IF ( psolver(1:9) == 'multigrid'  .OR.  psolver == 'sor' )  p = 0.0_wp
1501!
1502!--    Impose vortex with vertical axis on the initial velocity profile
1503       IF ( INDEX( initializing_actions, 'initialize_vortex' ) /= 0 )  THEN
1504          CALL init_rankine
1505       ENDIF
1506
1507!
1508!--    Impose temperature anomaly (advection test only) or warm air bubble
1509!--    close to surface
1510       IF ( INDEX( initializing_actions, 'initialize_ptanom' ) /= 0  .OR.  &
1511            INDEX( initializing_actions, 'initialize_bubble' ) /= 0  )  THEN
1512          CALL init_pt_anomaly
1513       ENDIF
1514       
1515!
1516!--    If required, change the surface temperature at the start of the 3D run
1517       IF ( pt_surface_initial_change /= 0.0_wp )  THEN
1518          pt(nzb,:,:) = pt(nzb,:,:) + pt_surface_initial_change
1519       ENDIF
1520
1521!
1522!--    If required, change the surface humidity/scalar at the start of the 3D
1523!--    run
1524       IF ( humidity  .AND.  q_surface_initial_change /= 0.0_wp )              &
1525          q(nzb,:,:) = q(nzb,:,:) + q_surface_initial_change
1526         
1527       IF ( passive_scalar .AND.  s_surface_initial_change /= 0.0_wp )         &
1528          s(nzb,:,:) = s(nzb,:,:) + s_surface_initial_change
1529       
1530
1531!
1532!--    Initialize old and new time levels.
1533       tpt_m = 0.0_wp; tu_m = 0.0_wp; tv_m = 0.0_wp; tw_m = 0.0_wp
1534       pt_p = pt; u_p = u; v_p = v; w_p = w
1535
1536       IF ( humidity  )  THEN
1537          tq_m = 0.0_wp
1538          q_p = q
1539       ENDIF
1540       
1541       IF ( passive_scalar )  THEN
1542          ts_m = 0.0_wp
1543          s_p  = s
1544       ENDIF       
1545
1546       IF ( debug_output )  CALL debug_message( 'initializing statistics, boundary conditions, etc.', 'end' )
1547
1548    ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data'  .OR.         &
1549             TRIM( initializing_actions ) == 'cyclic_fill' )                   &
1550    THEN
1551
1552       IF ( debug_output )  CALL debug_message( 'initializing in case of restart / cyclic_fill', 'start' )
1553!
1554!--    Initialize surface elements and its attributes, e.g. heat- and
1555!--    momentumfluxes, roughness, scaling parameters. As number of surface
1556!--    elements might be different between runs, e.g. in case of cyclic fill,
1557!--    and not all surface elements are read, surface elements need to be
1558!--    initialized before.
1559!--    Please note, in case of cyclic fill, surfaces should be initialized
1560!--    after restart data is read, else, individual settings of surface
1561!--    parameters will be overwritten from data of precursor run, hence,
1562!--    init_surfaces is called a second time after reading the restart data.
1563       CALL init_surfaces                       
1564!
1565!--    When reading data for cyclic fill of 3D prerun data files, read
1566!--    some of the global variables from the restart file which are required
1567!--    for initializing the inflow
1568       IF ( TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
1569
1570          DO  i = 0, io_blocks-1
1571             IF ( i == io_group )  THEN
1572                CALL rrd_read_parts_of_global
1573             ENDIF
1574#if defined( __parallel )
1575             CALL MPI_BARRIER( comm2d, ierr )
1576#endif
1577          ENDDO
1578
1579       ENDIF
1580
1581!
1582!--    Read processor specific binary data from restart file
1583       DO  i = 0, io_blocks-1
1584          IF ( i == io_group )  THEN
1585             CALL rrd_local
1586          ENDIF
1587#if defined( __parallel )
1588          CALL MPI_BARRIER( comm2d, ierr )
1589#endif
1590       ENDDO
1591!
1592!--    In case of cyclic fill, call init_surfaces a second time, so that
1593!--    surface properties such as heat fluxes are initialized as prescribed.
1594       IF ( TRIM( initializing_actions ) == 'cyclic_fill' )                    &
1595          CALL init_surfaces
1596
1597!
1598!--    In case of complex terrain and cyclic fill method as initialization,
1599!--    shift initial data in the vertical direction for each point in the
1600!--    x-y-plane depending on local surface height
1601       IF ( complex_terrain  .AND.                                             &
1602            TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
1603          DO  i = nxlg, nxrg
1604             DO  j = nysg, nyng
1605                nz_u_shift = get_topography_top_index_ji( j, i, 'u' )
1606                nz_v_shift = get_topography_top_index_ji( j, i, 'v' )
1607                nz_w_shift = get_topography_top_index_ji( j, i, 'w' )
1608                nz_s_shift = get_topography_top_index_ji( j, i, 's' )
1609
1610                u(nz_u_shift:nzt+1,j,i)  = u(0:nzt+1-nz_u_shift,j,i)               
1611
1612                v(nz_v_shift:nzt+1,j,i)  = v(0:nzt+1-nz_v_shift,j,i)
1613
1614                w(nz_w_shift:nzt+1,j,i)  = w(0:nzt+1-nz_w_shift,j,i)
1615
1616                p(nz_s_shift:nzt+1,j,i)  =  p(0:nzt+1-nz_s_shift,j,i)
1617                pt(nz_s_shift:nzt+1,j,i) = pt(0:nzt+1-nz_s_shift,j,i)
1618             ENDDO
1619          ENDDO
1620       ENDIF
1621
1622!
1623!--    Initialization of the turbulence recycling method
1624       IF ( TRIM( initializing_actions ) == 'cyclic_fill'  .AND.               &
1625            turbulent_inflow )  THEN
1626!
1627!--       First store the profiles to be used at the inflow.
1628!--       These profiles are the (temporally) and horizontally averaged vertical
1629!--       profiles from the prerun. Alternatively, prescribed profiles
1630!--       for u,v-components can be used.
1631          ALLOCATE( mean_inflow_profiles(nzb:nzt+1,1:num_mean_inflow_profiles) )
1632
1633          IF ( use_prescribed_profile_data )  THEN
1634             mean_inflow_profiles(:,1) = u_init            ! u
1635             mean_inflow_profiles(:,2) = v_init            ! v
1636          ELSE
1637             mean_inflow_profiles(:,1) = hom_sum(:,1,0)    ! u
1638             mean_inflow_profiles(:,2) = hom_sum(:,2,0)    ! v
1639          ENDIF
1640          mean_inflow_profiles(:,4) = hom_sum(:,4,0)       ! pt
1641          IF ( humidity )                                                      &
1642             mean_inflow_profiles(:,6) = hom_sum(:,41,0)   ! q
1643          IF ( passive_scalar )                                                &
1644             mean_inflow_profiles(:,7) = hom_sum(:,115,0)   ! s
1645!
1646!--       In case of complex terrain, determine vertical displacement at inflow
1647!--       boundary and adjust mean inflow profiles
1648          IF ( complex_terrain )  THEN
1649             IF ( nxlg <= 0 .AND. nxrg >= 0 .AND. nysg <= 0 .AND. nyng >= 0 )  THEN
1650                nz_u_shift_l = get_topography_top_index_ji( 0, 0, 'u' )
1651                nz_v_shift_l = get_topography_top_index_ji( 0, 0, 'v' )
1652                nz_w_shift_l = get_topography_top_index_ji( 0, 0, 'w' )
1653                nz_s_shift_l = get_topography_top_index_ji( 0, 0, 's' )
1654             ELSE
1655                nz_u_shift_l = 0
1656                nz_v_shift_l = 0
1657                nz_w_shift_l = 0
1658                nz_s_shift_l = 0
1659             ENDIF
1660
1661#if defined( __parallel )
1662             CALL MPI_ALLREDUCE(nz_u_shift_l, nz_u_shift, 1, MPI_INTEGER,      &
1663                                MPI_MAX, comm2d, ierr)
1664             CALL MPI_ALLREDUCE(nz_v_shift_l, nz_v_shift, 1, MPI_INTEGER,      &
1665                                MPI_MAX, comm2d, ierr)
1666             CALL MPI_ALLREDUCE(nz_w_shift_l, nz_w_shift, 1, MPI_INTEGER,      & 
1667                                MPI_MAX, comm2d, ierr)
1668             CALL MPI_ALLREDUCE(nz_s_shift_l, nz_s_shift, 1, MPI_INTEGER,      &
1669                                MPI_MAX, comm2d, ierr)
1670#else
1671             nz_u_shift = nz_u_shift_l
1672             nz_v_shift = nz_v_shift_l
1673             nz_w_shift = nz_w_shift_l
1674             nz_s_shift = nz_s_shift_l
1675#endif
1676
1677             mean_inflow_profiles(:,1) = 0.0_wp
1678             mean_inflow_profiles(nz_u_shift:nzt+1,1) = hom_sum(0:nzt+1-nz_u_shift,1,0)  ! u
1679
1680             mean_inflow_profiles(:,2) = 0.0_wp
1681             mean_inflow_profiles(nz_v_shift:nzt+1,2) = hom_sum(0:nzt+1-nz_v_shift,2,0)  ! v
1682
1683             mean_inflow_profiles(nz_s_shift:nzt+1,4) = hom_sum(0:nzt+1-nz_s_shift,4,0)  ! pt
1684
1685          ENDIF
1686
1687!
1688!--       If necessary, adjust the horizontal flow field to the prescribed
1689!--       profiles
1690          IF ( use_prescribed_profile_data )  THEN
1691             DO  i = nxlg, nxrg
1692                DO  j = nysg, nyng
1693                   DO  k = nzb, nzt+1
1694                      u(k,j,i) = u(k,j,i) - hom_sum(k,1,0) + u_init(k)
1695                      v(k,j,i) = v(k,j,i) - hom_sum(k,2,0) + v_init(k)
1696                   ENDDO
1697                ENDDO
1698             ENDDO
1699          ENDIF
1700
1701!
1702!--       Use these mean profiles at the inflow (provided that Dirichlet
1703!--       conditions are used)
1704          IF ( bc_dirichlet_l )  THEN
1705             DO  j = nysg, nyng
1706                DO  k = nzb, nzt+1
1707                   u(k,j,nxlg:-1)  = mean_inflow_profiles(k,1)
1708                   v(k,j,nxlg:-1)  = mean_inflow_profiles(k,2)
1709                   w(k,j,nxlg:-1)  = 0.0_wp
1710                   pt(k,j,nxlg:-1) = mean_inflow_profiles(k,4)
1711                   IF ( humidity )                                             &
1712                      q(k,j,nxlg:-1)  = mean_inflow_profiles(k,6)
1713                   IF ( passive_scalar )                                       &
1714                      s(k,j,nxlg:-1)  = mean_inflow_profiles(k,7)                     
1715                ENDDO
1716             ENDDO
1717          ENDIF
1718
1719!
1720!--       Calculate the damping factors to be used at the inflow. For a
1721!--       turbulent inflow the turbulent fluctuations have to be limited
1722!--       vertically because otherwise the turbulent inflow layer will grow
1723!--       in time.
1724          IF ( inflow_damping_height == 9999999.9_wp )  THEN
1725!
1726!--          Default: use the inversion height calculated by the prerun; if
1727!--          this is zero, inflow_damping_height must be explicitly
1728!--          specified.
1729             IF ( hom_sum(nzb+6,pr_palm,0) /= 0.0_wp )  THEN
1730                inflow_damping_height = hom_sum(nzb+6,pr_palm,0)
1731             ELSE
1732                WRITE( message_string, * ) 'inflow_damping_height must be ',   &
1733                     'explicitly specified because&the inversion height ',     &
1734                     'calculated by the prerun is zero.'
1735                CALL message( 'init_3d_model', 'PA0318', 1, 2, 0, 6, 0 )
1736             ENDIF
1737
1738          ENDIF
1739
1740          IF ( inflow_damping_width == 9999999.9_wp )  THEN
1741!
1742!--          Default for the transition range: one tenth of the undamped
1743!--          layer
1744             inflow_damping_width = 0.1_wp * inflow_damping_height
1745
1746          ENDIF
1747
1748          ALLOCATE( inflow_damping_factor(nzb:nzt+1) )
1749
1750          DO  k = nzb, nzt+1
1751
1752             IF ( zu(k) <= inflow_damping_height )  THEN
1753                inflow_damping_factor(k) = 1.0_wp
1754             ELSEIF ( zu(k) <= ( inflow_damping_height + inflow_damping_width ) )  THEN
1755                inflow_damping_factor(k) = 1.0_wp -                            &
1756                                           ( zu(k) - inflow_damping_height ) / &
1757                                           inflow_damping_width
1758             ELSE
1759                inflow_damping_factor(k) = 0.0_wp
1760             ENDIF
1761
1762          ENDDO
1763
1764       ENDIF
1765
1766!
1767!--    Inside buildings set velocities back to zero
1768       IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND.                &
1769            topography /= 'flat' )  THEN
1770!
1771!--       Inside buildings set velocities back to zero.
1772!--       Other scalars (pt, q, s, p, sa, ...) are ignored at present,
1773!--       maybe revise later.
1774          DO  i = nxlg, nxrg
1775             DO  j = nysg, nyng
1776                DO  k = nzb, nzt
1777                   u(k,j,i)     = MERGE( u(k,j,i), 0.0_wp,                     &
1778                                         BTEST( wall_flags_0(k,j,i), 1 ) )
1779                   v(k,j,i)     = MERGE( v(k,j,i), 0.0_wp,                     &
1780                                         BTEST( wall_flags_0(k,j,i), 2 ) )
1781                   w(k,j,i)     = MERGE( w(k,j,i), 0.0_wp,                     &
1782                                         BTEST( wall_flags_0(k,j,i), 3 ) )
1783                ENDDO
1784             ENDDO
1785          ENDDO
1786
1787       ENDIF
1788
1789!
1790!--    Calculate initial temperature field and other constants used in case
1791!--    of a sloping surface
1792       IF ( sloping_surface )  CALL init_slope
1793
1794!
1795!--    Initialize new time levels (only done in order to set boundary values
1796!--    including ghost points)
1797       pt_p = pt; u_p = u; v_p = v; w_p = w
1798       IF ( humidity )  THEN
1799          q_p = q
1800       ENDIF
1801       IF ( passive_scalar )  s_p  = s
1802!
1803!--    Allthough tendency arrays are set in prognostic_equations, they have
1804!--    have to be predefined here because they are used (but multiplied with 0)
1805!--    there before they are set.
1806       tpt_m = 0.0_wp; tu_m = 0.0_wp; tv_m = 0.0_wp; tw_m = 0.0_wp
1807       IF ( humidity )  THEN
1808          tq_m = 0.0_wp
1809       ENDIF
1810       IF ( passive_scalar )  ts_m  = 0.0_wp
1811
1812       IF ( debug_output )  CALL debug_message( 'initializing in case of restart / cyclic_fill', 'end' )
1813
1814    ELSE
1815!
1816!--    Actually this part of the programm should not be reached
1817       message_string = 'unknown initializing problem'
1818       CALL message( 'init_3d_model', 'PA0193', 1, 2, 0, 6, 0 )
1819    ENDIF
1820
1821
1822    IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1823!
1824!--    Initialize old timelevels needed for radiation boundary conditions
1825       IF ( bc_radiation_l )  THEN
1826          u_m_l(:,:,:) = u(:,:,1:2)
1827          v_m_l(:,:,:) = v(:,:,0:1)
1828          w_m_l(:,:,:) = w(:,:,0:1)
1829       ENDIF
1830       IF ( bc_radiation_r )  THEN
1831          u_m_r(:,:,:) = u(:,:,nx-1:nx)
1832          v_m_r(:,:,:) = v(:,:,nx-1:nx)
1833          w_m_r(:,:,:) = w(:,:,nx-1:nx)
1834       ENDIF
1835       IF ( bc_radiation_s )  THEN
1836          u_m_s(:,:,:) = u(:,0:1,:)
1837          v_m_s(:,:,:) = v(:,1:2,:)
1838          w_m_s(:,:,:) = w(:,0:1,:)
1839       ENDIF
1840       IF ( bc_radiation_n )  THEN
1841          u_m_n(:,:,:) = u(:,ny-1:ny,:)
1842          v_m_n(:,:,:) = v(:,ny-1:ny,:)
1843          w_m_n(:,:,:) = w(:,ny-1:ny,:)
1844       ENDIF
1845       
1846    ENDIF
1847
1848!
1849!-- Calculate the initial volume flow at the right and north boundary
1850    IF ( conserve_volume_flow )  THEN
1851
1852       IF ( use_prescribed_profile_data )  THEN
1853
1854          volume_flow_initial_l = 0.0_wp
1855          volume_flow_area_l    = 0.0_wp
1856
1857          IF ( nxr == nx )  THEN
1858             DO  j = nys, nyn
1859                DO  k = nzb+1, nzt
1860                   volume_flow_initial_l(1) = volume_flow_initial_l(1) +       &
1861                                              u_init(k) * dzw(k)               &
1862                                     * MERGE( 1.0_wp, 0.0_wp,                  &
1863                                              BTEST( wall_flags_0(k,j,nxr), 1 )&
1864                                            )
1865
1866                   volume_flow_area_l(1)    = volume_flow_area_l(1) + dzw(k)   &
1867                                     * MERGE( 1.0_wp, 0.0_wp,                  &
1868                                              BTEST( wall_flags_0(k,j,nxr), 1 )&
1869                                            )
1870                ENDDO
1871             ENDDO
1872          ENDIF
1873         
1874          IF ( nyn == ny )  THEN
1875             DO  i = nxl, nxr
1876                DO  k = nzb+1, nzt
1877                   volume_flow_initial_l(2) = volume_flow_initial_l(2) +       &
1878                                              v_init(k) * dzw(k)               &       
1879                                     * MERGE( 1.0_wp, 0.0_wp,                  &
1880                                              BTEST( wall_flags_0(k,nyn,i), 2 )&
1881                                            )
1882                   volume_flow_area_l(2)    = volume_flow_area_l(2) + dzw(k)   &       
1883                                     * MERGE( 1.0_wp, 0.0_wp,                  &
1884                                              BTEST( wall_flags_0(k,nyn,i), 2 )&
1885                                            )
1886                ENDDO
1887             ENDDO
1888          ENDIF
1889
1890#if defined( __parallel )
1891          CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&
1892                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
1893          CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1),      &
1894                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
1895
1896#else
1897          volume_flow_initial = volume_flow_initial_l
1898          volume_flow_area    = volume_flow_area_l
1899#endif 
1900
1901       ELSEIF ( TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
1902
1903          volume_flow_initial_l = 0.0_wp
1904          volume_flow_area_l    = 0.0_wp
1905
1906          IF ( nxr == nx )  THEN
1907             DO  j = nys, nyn
1908                DO  k = nzb+1, nzt
1909                   volume_flow_initial_l(1) = volume_flow_initial_l(1) +       &
1910                                              hom_sum(k,1,0) * dzw(k)          &
1911                                     * MERGE( 1.0_wp, 0.0_wp,                  &
1912                                              BTEST( wall_flags_0(k,j,nx), 1 ) &
1913                                            )
1914                   volume_flow_area_l(1)    = volume_flow_area_l(1) + dzw(k)   &
1915                                     * MERGE( 1.0_wp, 0.0_wp,                  &
1916                                              BTEST( wall_flags_0(k,j,nx), 1 ) &
1917                                            )
1918                ENDDO
1919             ENDDO
1920          ENDIF
1921         
1922          IF ( nyn == ny )  THEN
1923             DO  i = nxl, nxr
1924                DO  k = nzb+1, nzt
1925                   volume_flow_initial_l(2) = volume_flow_initial_l(2) +       &
1926                                              hom_sum(k,2,0) * dzw(k)          &       
1927                                     * MERGE( 1.0_wp, 0.0_wp,                  &
1928                                              BTEST( wall_flags_0(k,ny,i), 2 ) &
1929                                            )
1930                   volume_flow_area_l(2)    = volume_flow_area_l(2) + dzw(k)   &       
1931                                     * MERGE( 1.0_wp, 0.0_wp,                  &
1932                                              BTEST( wall_flags_0(k,ny,i), 2 ) &
1933                                            )
1934                ENDDO
1935             ENDDO
1936          ENDIF
1937
1938#if defined( __parallel )
1939          CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&
1940                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
1941          CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1),      &
1942                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
1943
1944#else
1945          volume_flow_initial = volume_flow_initial_l
1946          volume_flow_area    = volume_flow_area_l
1947#endif 
1948
1949       ELSEIF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1950
1951          volume_flow_initial_l = 0.0_wp
1952          volume_flow_area_l    = 0.0_wp
1953
1954          IF ( nxr == nx )  THEN
1955             DO  j = nys, nyn
1956                DO  k = nzb+1, nzt
1957                   volume_flow_initial_l(1) = volume_flow_initial_l(1) +       &
1958                                              u(k,j,nx) * dzw(k)               &
1959                                     * MERGE( 1.0_wp, 0.0_wp,                  &
1960                                              BTEST( wall_flags_0(k,j,nx), 1 ) &
1961                                            )
1962                   volume_flow_area_l(1)    = volume_flow_area_l(1) + dzw(k)   &
1963                                     * MERGE( 1.0_wp, 0.0_wp,                  &
1964                                              BTEST( wall_flags_0(k,j,nx), 1 ) &
1965                                            )
1966                ENDDO
1967             ENDDO
1968          ENDIF
1969         
1970          IF ( nyn == ny )  THEN
1971             DO  i = nxl, nxr
1972                DO  k = nzb+1, nzt
1973                   volume_flow_initial_l(2) = volume_flow_initial_l(2) +       &
1974                                              v(k,ny,i) * dzw(k)               &       
1975                                     * MERGE( 1.0_wp, 0.0_wp,                  &
1976                                              BTEST( wall_flags_0(k,ny,i), 2 ) &
1977                                            )
1978                   volume_flow_area_l(2)    = volume_flow_area_l(2) + dzw(k)   &       
1979                                     * MERGE( 1.0_wp, 0.0_wp,                  &
1980                                              BTEST( wall_flags_0(k,ny,i), 2 ) &
1981                                            )
1982                ENDDO
1983             ENDDO
1984          ENDIF
1985
1986#if defined( __parallel )
1987          CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&
1988                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
1989          CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1),      &
1990                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
1991
1992#else
1993          volume_flow_initial = volume_flow_initial_l
1994          volume_flow_area    = volume_flow_area_l
1995#endif 
1996
1997       ENDIF
1998
1999!
2000!--    In case of 'bulk_velocity' mode, volume_flow_initial is calculated
2001!--    from u|v_bulk instead
2002       IF ( TRIM( conserve_volume_flow_mode ) == 'bulk_velocity' )  THEN
2003          volume_flow_initial(1) = u_bulk * volume_flow_area(1)
2004          volume_flow_initial(2) = v_bulk * volume_flow_area(2)
2005       ENDIF
2006
2007    ENDIF
2008!
2009!-- Finally, if random_heatflux is set, disturb shf at horizontal
2010!-- surfaces. Actually, this should be done in surface_mod, where all other
2011!-- initializations of surface quantities are done. However, this
2012!-- would create a ring dependency, hence, it is done here. Maybe delete
2013!-- disturb_heatflux and tranfer the respective code directly into the
2014!-- initialization in surface_mod.         
2015    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
2016         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
2017 
2018       IF ( use_surface_fluxes  .AND.  constant_heatflux  .AND.                &
2019            random_heatflux )  THEN
2020          IF ( surf_def_h(0)%ns >= 1 )  CALL disturb_heatflux( surf_def_h(0) )
2021          IF ( surf_lsm_h%ns    >= 1 )  CALL disturb_heatflux( surf_lsm_h    )
2022          IF ( surf_usm_h%ns    >= 1 )  CALL disturb_heatflux( surf_usm_h    )
2023       ENDIF
2024    ENDIF
2025
2026!
2027!-- Compute total sum of grid points and the mean surface level height for each
2028!-- statistic region. These are mainly used for horizontal averaging of
2029!-- turbulence statistics.
2030!-- ngp_2dh: number of grid points of a horizontal cross section through the
2031!--          respective statistic region
2032!-- ngp_3d:  number of grid points of the respective statistic region
2033    ngp_2dh_outer_l   = 0
2034    ngp_2dh_outer     = 0
2035    ngp_2dh_s_inner_l = 0
2036    ngp_2dh_s_inner   = 0
2037    ngp_2dh_l         = 0
2038    ngp_2dh           = 0
2039    ngp_3d_inner_l    = 0.0_wp
2040    ngp_3d_inner      = 0
2041    ngp_3d            = 0
2042    ngp_sums          = ( nz + 2 ) * ( pr_palm + max_pr_user )
2043
2044    mean_surface_level_height   = 0.0_wp
2045    mean_surface_level_height_l = 0.0_wp
2046!
2047!-- To do: New concept for these non-topography grid points!
2048    DO  sr = 0, statistic_regions
2049       DO  i = nxl, nxr
2050          DO  j = nys, nyn
2051             IF ( rmask(j,i,sr) == 1.0_wp )  THEN
2052!
2053!--             All xy-grid points
2054                ngp_2dh_l(sr) = ngp_2dh_l(sr) + 1
2055!
2056!--             Determine mean surface-level height. In case of downward-
2057!--             facing walls are present, more than one surface level exist.
2058!--             In this case, use the lowest surface-level height.
2059                IF ( surf_def_h(0)%start_index(j,i) <=                         &
2060                     surf_def_h(0)%end_index(j,i) )  THEN
2061                   m = surf_def_h(0)%start_index(j,i)
2062                   k = surf_def_h(0)%k(m)
2063                   mean_surface_level_height_l(sr) =                           &
2064                                       mean_surface_level_height_l(sr) + zw(k-1)
2065                ENDIF
2066                IF ( surf_lsm_h%start_index(j,i) <=                            &
2067                     surf_lsm_h%end_index(j,i) )  THEN
2068                   m = surf_lsm_h%start_index(j,i)
2069                   k = surf_lsm_h%k(m)
2070                   mean_surface_level_height_l(sr) =                           &
2071                                       mean_surface_level_height_l(sr) + zw(k-1)
2072                ENDIF
2073                IF ( surf_usm_h%start_index(j,i) <=                            &
2074                     surf_usm_h%end_index(j,i) )  THEN
2075                   m = surf_usm_h%start_index(j,i)
2076                   k = surf_usm_h%k(m)
2077                   mean_surface_level_height_l(sr) =                           &
2078                                       mean_surface_level_height_l(sr) + zw(k-1)
2079                ENDIF
2080
2081                k_surf = k - 1
2082
2083                DO  k = nzb, nzt+1
2084!
2085!--                xy-grid points above topography
2086                   ngp_2dh_outer_l(k,sr) = ngp_2dh_outer_l(k,sr)     +         &
2087                                  MERGE( 1, 0, BTEST( wall_flags_0(k,j,i), 24 ) )
2088
2089                   ngp_2dh_s_inner_l(k,sr) = ngp_2dh_s_inner_l(k,sr) +         &
2090                                  MERGE( 1, 0, BTEST( wall_flags_0(k,j,i), 22 ) )
2091
2092                ENDDO
2093!
2094!--             All grid points of the total domain above topography
2095                ngp_3d_inner_l(sr) = ngp_3d_inner_l(sr) + ( nz - k_surf + 2 )
2096
2097
2098
2099             ENDIF
2100          ENDDO
2101       ENDDO
2102    ENDDO
2103
2104    sr = statistic_regions + 1
2105#if defined( __parallel )
2106    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2107    CALL MPI_ALLREDUCE( ngp_2dh_l(0), ngp_2dh(0), sr, MPI_INTEGER, MPI_SUM,    &
2108                        comm2d, ierr )
2109    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2110    CALL MPI_ALLREDUCE( ngp_2dh_outer_l(0,0), ngp_2dh_outer(0,0), (nz+2)*sr,   &
2111                        MPI_INTEGER, MPI_SUM, comm2d, ierr )
2112    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2113    CALL MPI_ALLREDUCE( ngp_2dh_s_inner_l(0,0), ngp_2dh_s_inner(0,0),          &
2114                        (nz+2)*sr, MPI_INTEGER, MPI_SUM, comm2d, ierr )
2115    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2116    CALL MPI_ALLREDUCE( ngp_3d_inner_l(0), ngp_3d_inner_tmp(0), sr, MPI_REAL,  &
2117                        MPI_SUM, comm2d, ierr )
2118    ngp_3d_inner = INT( ngp_3d_inner_tmp, KIND = SELECTED_INT_KIND( 18 ) )
2119    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2120    CALL MPI_ALLREDUCE( mean_surface_level_height_l(0),                        &
2121                        mean_surface_level_height(0), sr, MPI_REAL,            &
2122                        MPI_SUM, comm2d, ierr )
2123    mean_surface_level_height = mean_surface_level_height / REAL( ngp_2dh )
2124#else
2125    ngp_2dh         = ngp_2dh_l
2126    ngp_2dh_outer   = ngp_2dh_outer_l
2127    ngp_2dh_s_inner = ngp_2dh_s_inner_l
2128    ngp_3d_inner    = INT( ngp_3d_inner_l, KIND = SELECTED_INT_KIND( 18 ) )
2129    mean_surface_level_height = mean_surface_level_height_l / REAL( ngp_2dh_l )
2130#endif
2131
2132    ngp_3d = INT ( ngp_2dh, KIND = SELECTED_INT_KIND( 18 ) ) * &
2133             INT ( (nz + 2 ), KIND = SELECTED_INT_KIND( 18 ) )
2134
2135!
2136!-- Set a lower limit of 1 in order to avoid zero divisions in flow_statistics,
2137!-- buoyancy, etc. A zero value will occur for cases where all grid points of
2138!-- the respective subdomain lie below the surface topography
2139    ngp_2dh_outer   = MAX( 1, ngp_2dh_outer(:,:)   ) 
2140    ngp_3d_inner    = MAX( INT(1, KIND = SELECTED_INT_KIND( 18 )),             &
2141                           ngp_3d_inner(:) )
2142    ngp_2dh_s_inner = MAX( 1, ngp_2dh_s_inner(:,:) ) 
2143
2144    DEALLOCATE( mean_surface_level_height_l, ngp_2dh_l, ngp_2dh_outer_l,       &
2145                ngp_3d_inner_l, ngp_3d_inner_tmp )
2146!
2147!-- Initialize surface forcing corresponding to large-scale forcing. Therein,
2148!-- initialize heat-fluxes, etc. via datatype. Revise it later!
2149    IF ( large_scale_forcing .AND. lsf_surf )  THEN
2150       IF ( use_surface_fluxes  .AND.  constant_heatflux )  THEN
2151          CALL ls_forcing_surf ( simulated_time )
2152       ENDIF
2153    ENDIF
2154!
2155!-- Initializae 3D offline nesting in COSMO model and read data from
2156!-- external NetCDF file.
2157    IF ( nesting_offline )  CALL nesting_offl_init
2158!
2159!-- Initialize quantities for special advections schemes
2160    CALL init_advec
2161
2162!
2163!-- Impose random perturbation on the horizontal velocity field and then
2164!-- remove the divergences from the velocity field at the initial stage
2165    IF ( create_disturbances  .AND.  disturbance_energy_limit /= 0.0_wp  .AND. &
2166         TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
2167         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
2168
2169       IF ( debug_output )  CALL debug_message( 'creating disturbances + applying pressure solver', 'start' )
2170!
2171!--    Needed for both disturb_field and pres
2172!$ACC DATA &
2173!$ACC CREATE(tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
2174!$ACC COPY(u(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
2175!$ACC COPY(v(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
2176
2177       CALL disturb_field( 'u', tend, u )
2178       CALL disturb_field( 'v', tend, v )
2179
2180!$ACC DATA &
2181!$ACC CREATE(d(nzb+1:nzt,nys:nyn,nxl:nxr)) &
2182!$ACC COPY(w(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
2183!$ACC COPY(p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
2184!$ACC COPYIN(rho_air(nzb:nzt+1), rho_air_zw(nzb:nzt+1)) &
2185!$ACC COPYIN(ddzu(1:nzt+1), ddzw(1:nzt+1)) &
2186!$ACC COPYIN(wall_flags_0(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
2187!$ACC COPYIN(bc_h(0:1)) &
2188!$ACC COPYIN(bc_h(0)%i(1:bc_h(0)%ns)) &
2189!$ACC COPYIN(bc_h(0)%j(1:bc_h(0)%ns)) &
2190!$ACC COPYIN(bc_h(0)%k(1:bc_h(0)%ns)) &
2191!$ACC COPYIN(bc_h(1)%i(1:bc_h(1)%ns)) &
2192!$ACC COPYIN(bc_h(1)%j(1:bc_h(1)%ns)) &
2193!$ACC COPYIN(bc_h(1)%k(1:bc_h(1)%ns))
2194
2195       n_sor = nsor_ini
2196       CALL pres
2197       n_sor = nsor
2198
2199!$ACC END DATA
2200!$ACC END DATA
2201
2202       IF ( debug_output )  CALL debug_message( 'creating disturbances + applying pressure solver', 'end' )
2203
2204    ENDIF
2205
2206    IF ( .NOT. ocean_mode )  THEN
2207
2208       ALLOCATE( hyp(nzb:nzt+1) )
2209       ALLOCATE( d_exner(nzb:nzt+1) )
2210       ALLOCATE( exner(nzb:nzt+1) )
2211       ALLOCATE( hyrho(nzb:nzt+1) )
2212!
2213!--    Check temperature in case of too large domain height
2214       DO  k = nzb, nzt+1
2215          IF ( ( pt_surface * exner_function(surface_pressure * 100.0_wp) - g/c_p * zu(k) ) < 0.0_wp )  THEN
2216             WRITE( message_string, * )  'absolute temperature < 0.0 at zu(', k, &
2217                                         ') = ', zu(k)
2218             CALL message( 'init_3d_model', 'PA0142', 1, 2, 0, 6, 0 )
2219          ENDIF
2220       ENDDO
2221
2222!
2223!--    Calculate vertical profile of the hydrostatic pressure (hyp)
2224       hyp    = barometric_formula(zu, pt_surface * exner_function(surface_pressure * 100.0_wp), surface_pressure * 100.0_wp)
2225       d_exner = exner_function_invers(hyp)
2226       exner = 1.0_wp / exner_function_invers(hyp)
2227       hyrho  = ideal_gas_law_rho_pt(hyp, pt_init)
2228!
2229!--    Compute reference density
2230       rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp, pt_surface * exner_function(surface_pressure * 100.0_wp))
2231
2232    ENDIF
2233
2234!
2235!-- If required, initialize particles
2236    IF ( agents_active )  CALL mas_init
2237!
2238!-- Initialization of synthetic turbulence generator
2239    IF ( use_syn_turb_gen )  CALL stg_init
2240!
2241!-- Initializing actions for all other modules
2242    CALL module_interface_init
2243!
2244!-- Initialize surface layer (done after LSM as roughness length are required
2245!-- for initialization
2246    IF ( constant_flux_layer )  CALL init_surface_layer_fluxes
2247!
2248!-- Initialize surface data output
2249    IF ( surface_output )  CALL surface_data_output_init
2250!
2251!-- Initialize the ws-scheme.   
2252    IF ( ws_scheme_sca .OR. ws_scheme_mom )  CALL ws_init
2253!
2254!-- Perform post-initializing checks for all other modules
2255    CALL module_interface_init_checks
2256
2257!
2258!-- Setting weighting factors for calculation of perturbation pressure
2259!-- and turbulent quantities from the RK substeps
2260    IF ( TRIM(timestep_scheme) == 'runge-kutta-3' )  THEN      ! for RK3-method
2261
2262       weight_substep(1) = 1._wp/6._wp
2263       weight_substep(2) = 3._wp/10._wp
2264       weight_substep(3) = 8._wp/15._wp
2265
2266       weight_pres(1)    = 1._wp/3._wp
2267       weight_pres(2)    = 5._wp/12._wp
2268       weight_pres(3)    = 1._wp/4._wp
2269
2270    ELSEIF ( TRIM(timestep_scheme) == 'runge-kutta-2' )  THEN  ! for RK2-method
2271
2272       weight_substep(1) = 1._wp/2._wp
2273       weight_substep(2) = 1._wp/2._wp
2274         
2275       weight_pres(1)    = 1._wp/2._wp
2276       weight_pres(2)    = 1._wp/2._wp       
2277
2278    ELSE                                     ! for Euler-method
2279
2280       weight_substep(1) = 1.0_wp     
2281       weight_pres(1)    = 1.0_wp                   
2282
2283    ENDIF
2284
2285!
2286!-- Initialize Rayleigh damping factors
2287    rdf    = 0.0_wp
2288    rdf_sc = 0.0_wp
2289    IF ( rayleigh_damping_factor /= 0.0_wp )  THEN
2290
2291       IF (  .NOT.  ocean_mode )  THEN
2292          DO  k = nzb+1, nzt
2293             IF ( zu(k) >= rayleigh_damping_height )  THEN
2294                rdf(k) = rayleigh_damping_factor *                             &
2295                      ( SIN( pi * 0.5_wp * ( zu(k) - rayleigh_damping_height ) &
2296                             / ( zu(nzt) - rayleigh_damping_height ) )         &
2297                      )**2
2298             ENDIF
2299          ENDDO
2300       ELSE
2301!
2302!--       In ocean mode, rayleigh damping is applied in the lower part of the
2303!--       model domain
2304          DO  k = nzt, nzb+1, -1
2305             IF ( zu(k) <= rayleigh_damping_height )  THEN
2306                rdf(k) = rayleigh_damping_factor *                             &
2307                      ( SIN( pi * 0.5_wp * ( rayleigh_damping_height - zu(k) ) &
2308                             / ( rayleigh_damping_height - zu(nzb+1) ) )       &
2309                      )**2
2310             ENDIF
2311          ENDDO
2312       ENDIF
2313
2314    ENDIF
2315    IF ( scalar_rayleigh_damping )  rdf_sc = rdf
2316
2317!
2318!-- Initialize the starting level and the vertical smoothing factor used for
2319!-- the external pressure gradient
2320    dp_smooth_factor = 1.0_wp
2321    IF ( dp_external )  THEN
2322!
2323!--    Set the starting level dp_level_ind_b only if it has not been set before
2324!--    (e.g. in init_grid).
2325       IF ( dp_level_ind_b == 0 )  THEN
2326          ind_array = MINLOC( ABS( dp_level_b - zu ) )
2327          dp_level_ind_b = ind_array(1) - 1 + nzb
2328                                        ! MINLOC uses lower array bound 1
2329       ENDIF
2330       IF ( dp_smooth )  THEN
2331          dp_smooth_factor(:dp_level_ind_b) = 0.0_wp
2332          DO  k = dp_level_ind_b+1, nzt
2333             dp_smooth_factor(k) = 0.5_wp * ( 1.0_wp + SIN( pi *               &
2334                        ( REAL( k - dp_level_ind_b, KIND=wp ) /                &
2335                          REAL( nzt - dp_level_ind_b, KIND=wp ) - 0.5_wp ) ) )
2336          ENDDO
2337       ENDIF
2338    ENDIF
2339
2340!
2341!-- Initialize damping zone for the potential temperature in case of
2342!-- non-cyclic lateral boundaries. The damping zone has the maximum value
2343!-- at the inflow boundary and decreases to zero at pt_damping_width.
2344    ptdf_x = 0.0_wp
2345    ptdf_y = 0.0_wp
2346    IF ( bc_lr_dirrad )  THEN
2347       DO  i = nxl, nxr
2348          IF ( ( i * dx ) < pt_damping_width )  THEN
2349             ptdf_x(i) = pt_damping_factor * ( SIN( pi * 0.5_wp *              &
2350                            REAL( pt_damping_width - i * dx, KIND=wp ) / (     &
2351                            REAL( pt_damping_width, KIND=wp ) ) ) )**2 
2352          ENDIF
2353       ENDDO
2354    ELSEIF ( bc_lr_raddir )  THEN
2355       DO  i = nxl, nxr
2356          IF ( ( i * dx ) > ( nx * dx - pt_damping_width ) )  THEN
2357             ptdf_x(i) = pt_damping_factor *                                   &
2358                         SIN( pi * 0.5_wp *                                    &
2359                                 ( ( i - nx ) * dx + pt_damping_width ) /      &
2360                                 REAL( pt_damping_width, KIND=wp ) )**2
2361          ENDIF
2362       ENDDO
2363    ELSEIF ( bc_ns_dirrad )  THEN
2364       DO  j = nys, nyn
2365          IF ( ( j * dy ) > ( ny * dy - pt_damping_width ) )  THEN
2366             ptdf_y(j) = pt_damping_factor *                                   &
2367                         SIN( pi * 0.5_wp *                                    &
2368                                 ( ( j - ny ) * dy + pt_damping_width ) /      &
2369                                 REAL( pt_damping_width, KIND=wp ) )**2
2370          ENDIF
2371       ENDDO
2372    ELSEIF ( bc_ns_raddir )  THEN
2373       DO  j = nys, nyn
2374          IF ( ( j * dy ) < pt_damping_width )  THEN
2375             ptdf_y(j) = pt_damping_factor *                                   &
2376                         SIN( pi * 0.5_wp *                                    &
2377                                ( pt_damping_width - j * dy ) /                &
2378                                REAL( pt_damping_width, KIND=wp ) )**2
2379          ENDIF
2380       ENDDO
2381    ENDIF
2382
2383!
2384!-- Input binary data file is not needed anymore. This line must be placed
2385!-- after call of user_init!
2386    CALL close_file( 13 )
2387!
2388!-- In case of nesting, put an barrier to assure that all parent and child
2389!-- domains finished initialization.
2390#if defined( __parallel )
2391    IF ( nested_run )  CALL MPI_BARRIER( MPI_COMM_WORLD, ierr )
2392#endif
2393
2394
2395    CALL location_message( 'model initialization', 'finished' )
2396
2397 END SUBROUTINE init_3d_model
Note: See TracBrowser for help on using the repository browser.