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

Last change on this file since 2233 was 2233, checked in by suehring, 4 years ago

last commit documented

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