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

Last change on this file since 1738 was 1738, checked in by raasch, 6 years ago

bugfixes for calculations in statistical regions which do not contain grid points in the lowest vertical levels, mean surface level height considered in the calculation of the characteristic vertical velocity

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