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

Last change on this file since 1484 was 1484, checked in by kanani, 9 years ago

New:
---
Subroutine init_plant_canopy added to module plant_canopy_model_mod. (plant_canopy_model)
Alternative method for lad-profile construction added, also, new parameters added.
(header, package_parin, plant_canopy_model, read_var_list, write_var_list)
plant_canopy_model-dependency added to several subroutines. (Makefile)
New package/namelist canopy_par for canopy-related parameters added. (package_parin)

Changed:
---
Code structure of the plant canopy model changed, all canopy-model related code
combined to module plant_canopy_model_mod. (check_parameters, init_3d_model,
modules, timestep)
Module plant_canopy_model_mod added in USE-lists of some subroutines. (check_parameters,
header, init_3d_model, package_parin, read_var_list, user_init_plant_canopy, write_var_list)
Canopy initialization moved to new subroutine init_plant_canopy. (check_parameters,
init_3d_model, plant_canopy_model)
Calculation of canopy timestep-criterion removed, instead, the canopy
drag is now directly limited in the calculation of the canopy tendency terms.
(plant_canopy_model, timestep)
Some parameters renamed. (check_parameters, header, init_plant_canopy,
plant_canopy_model, read_var_list, write_var_list)
Unnecessary 3d-arrays removed. (init_plant_canopy, plant_canopy_model, user_init_plant_canopy)
Parameter checks regarding canopy initialization added. (check_parameters)
All canopy steering parameters moved from namelist inipar to canopy_par. (package_parin, parin)
Some redundant MPI communication removed. (init_plant_canopy)

Bugfix:
---
Missing KIND-attribute for REAL constant added. (check_parameters)
DO-WHILE-loop for lad-profile output restricted. (header)
Removed double-listing of use_upstream_for_tke in ONLY-list of module
control_parameters. (prognostic_equations)

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