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

Last change on this file since 4185 was 4185, checked in by oliver.maas, 2 years ago

For initializing_actions = ' cyclic_fill': Overwrite u_init, v_init, pt_init, q_init and s_init with the (temporally) and horizontally averaged vertical profiles from the end of the prerun, because these profiles shall be used as the basic state for the rayleigh damping and the pt_damping.

  • Property svn:keywords set to Id
  • Property svn:mergeinfo set to False
    /palm/branches/chemistry/SOURCE/init_3d_model.f902047-3190,​3218-3297
    /palm/branches/forwind/SOURCE/init_3d_model.f901564-1913
    /palm/branches/mosaik_M2/init_3d_model.f902360-3471
    /palm/branches/palm4u/SOURCE/init_3d_model.f902540-2692
    /palm/branches/rans/SOURCE/init_3d_model.f902078-3128
    /palm/branches/resler/SOURCE/init_3d_model.f902023-3605
    /palm/branches/salsa/SOURCE/init_3d_model.f902503-3581
File size: 80.7 KB
Line 
1!> @file init_3d_model.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2019 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: init_3d_model.f90 4185 2019-08-23 13:49:38Z oliver.maas $
27! For initializing_actions = ' cyclic_fill':
28! Overwrite u_init, v_init, pt_init, q_init and s_init with the
29! (temporally) and horizontally averaged vertical profiles from the end
30! of the prerun, because these profiles shall be used as the basic state
31! for the rayleigh damping and the pt_damping.
32!
33! 4182 2019-08-22 15:20:23Z scharf
34! Corrected "Former revisions" section
35!
36! 4168 2019-08-16 13:50:17Z suehring
37! Replace function get_topography_top_index by topo_top_ind
38!
39! 4151 2019-08-09 08:24:30Z suehring
40! Add netcdf directive around input calls (fix for last commit)
41!
42! 4150 2019-08-08 20:00:47Z suehring
43! Input of additional surface variables independent on land- or urban-surface
44! model
45!
46! 4131 2019-08-02 11:06:18Z monakurppa
47! Allocate sums and sums_l to allow profile output for salsa variables.
48!
49! 4130 2019-08-01 13:04:13Z suehring
50! Effectively reduce 3D initialization to 1D initial profiles. This is because
51! 3D initialization produces structures in the w-component that are correlated
52! with the processor grid for some unknown reason 
53!
54! 4090 2019-07-11 15:06:47Z Giersch
55! Unused variables removed
56!
57! 4088 2019-07-11 13:57:56Z Giersch
58! Pressure and density profile calculations revised using basic functions
59!
60! 4048 2019-06-21 21:00:21Z knoop
61! Further modularization of particle code components
62!
63! 4017 2019-06-06 12:16:46Z schwenkel
64! Convert most location messages to debug messages to reduce output in
65! job logfile to a minimum
66!
67!
68! unused variable removed
69!
70! 3937 2019-04-29 15:09:07Z suehring
71! Move initialization of synthetic turbulence generator behind initialization
72! of offline nesting. Remove call for stg_adjust, as this is now already done
73! in stg_init.
74!
75! 3900 2019-04-16 15:17:43Z suehring
76! Fix problem with LOD = 2 initialization
77!
78! 3885 2019-04-11 11:29:34Z kanani
79! Changes related to global restructuring of location messages and introduction
80! of additional debug messages
81!
82! 3849 2019-04-01 16:35:16Z knoop
83! Move initialization of rmask before initializing user_init_arrays
84!
85! 3711 2019-01-31 13:44:26Z knoop
86! Introduced module_interface_init_checks for post-init checks in modules
87!
88! 3700 2019-01-26 17:03:42Z knoop
89! Some interface calls moved to module_interface + cleanup
90!
91! 3648 2019-01-02 16:35:46Z suehring
92! Rename subroutines for surface-data output
93!
94! Revision 1.1  1998/03/09 16:22:22  raasch
95! Initial revision
96!
97!
98! Description:
99! ------------
100!> Allocation of arrays and initialization of the 3D model via
101!> a) pre-run the 1D model
102!> or
103!> b) pre-set constant linear profiles
104!> or
105!> c) read values of a previous run
106!------------------------------------------------------------------------------!
107 SUBROUTINE init_3d_model
108
109
110    USE advec_ws
111
112    USE arrays_3d
113
114    USE basic_constants_and_equations_mod,                                     &
115        ONLY:  c_p, g, l_v, pi, exner_function, exner_function_invers,         &
116               ideal_gas_law_rho, ideal_gas_law_rho_pt, barometric_formula
117
118    USE bulk_cloud_model_mod,                                                  &
119        ONLY:  bulk_cloud_model
120
121    USE chem_modules,                                                          &
122        ONLY:  max_pr_cs ! ToDo: this dependency needs to be removed cause it is ugly #new_dom
123
124    USE control_parameters
125
126    USE grid_variables,                                                        &
127        ONLY:  dx, dy, ddx2_mg, ddy2_mg
128
129    USE indices
130
131    USE kinds
132 
133    USE lsf_nudging_mod,                                                       &
134        ONLY:  ls_forcing_surf
135
136    USE model_1d_mod,                                                          &
137        ONLY:  init_1d_model, l1d, u1d, v1d
138
139    USE module_interface,                                                      &
140        ONLY:  module_interface_init_arrays,                                   &
141               module_interface_init,                                          &
142               module_interface_init_checks
143
144    USE multi_agent_system_mod,                                                &
145        ONLY:  agents_active, mas_init
146
147    USE netcdf_interface,                                                      &
148        ONLY:  dots_max
149
150    USE netcdf_data_input_mod,                                                 &
151        ONLY:  char_fill,                                                      &
152               check_existence,                                                &
153               close_input_file,                                               &
154               get_attribute,                                                  &
155               get_variable,                                                   &
156               init_3d,                                                        &
157               input_pids_static,                                              &
158               inquire_num_variables,                                          &
159               inquire_variable_names,                                         &
160               input_file_static,                                              &
161               netcdf_data_input_init_3d,                                      &
162               open_read_file,                                                 &
163               real_2d
164               
165    USE nesting_offl_mod,                                                      &
166        ONLY:  nesting_offl_init
167
168    USE pegrid
169
170#if defined( __parallel )
171    USE pmc_interface,                                                         &
172        ONLY:  nested_run
173#endif
174
175    USE random_function_mod
176
177    USE random_generator_parallel,                                             &
178        ONLY:  init_parallel_random_generator
179
180    USE read_restart_data_mod,                                                 &
181        ONLY:  rrd_read_parts_of_global, rrd_local
182
183    USE statistics,                                                            &
184        ONLY:  hom, hom_sum, mean_surface_level_height, pr_palm, rmask,        &
185               statistic_regions, sums, sums_divnew_l, sums_divold_l, sums_l,  &
186               sums_l_l, sums_wsts_bc_l, ts_value,                             &
187               weight_pres, weight_substep
188
189    USE synthetic_turbulence_generator_mod,                                    &
190        ONLY:  stg_init, use_syn_turb_gen
191
192    USE surface_layer_fluxes_mod,                                              &
193        ONLY:  init_surface_layer_fluxes
194
195    USE surface_mod,                                                           &
196        ONLY :  init_single_surface_properties,                                &
197                init_surface_arrays,                                           &
198                init_surfaces,                                                 &
199                surf_def_h,                                                    &
200                surf_def_v,                                                    &
201                surf_lsm_h,                                                    &
202                surf_usm_h
203
204#if defined( _OPENACC )
205    USE surface_mod,                                                           &
206        ONLY :  bc_h
207#endif
208
209    USE surface_data_output_mod,                                               &
210        ONLY:  surface_data_output_init
211
212    USE transpose_indices
213
214    IMPLICIT NONE
215
216    CHARACTER(LEN=50), DIMENSION(:), ALLOCATABLE ::  vars_pids_static  !< variable names in static input file
217   
218    INTEGER(iwp) ::  i                    !< grid index in x direction
219    INTEGER(iwp) ::  ind_array(1)         !< dummy used to determine start index for external pressure forcing
220    INTEGER(iwp) ::  j                    !< grid index in y direction
221    INTEGER(iwp) ::  k                    !< grid index in z direction
222    INTEGER(iwp) ::  k_surf               !< surface level index
223    INTEGER(iwp) ::  l                    !< running index over surface orientation   
224    INTEGER(iwp) ::  m                    !< index of surface element in surface data type   
225    INTEGER(iwp) ::  num_vars_pids_static !< number of variables in static input file
226    INTEGER(iwp) ::  nz_u_shift           !< topography-top index on u-grid, used to vertically shift initial profiles
227    INTEGER(iwp) ::  nz_v_shift           !< topography-top index on v-grid, used to vertically shift initial profiles
228    INTEGER(iwp) ::  nz_w_shift           !< topography-top index on w-grid, used to vertically shift initial profiles
229    INTEGER(iwp) ::  nz_s_shift           !< topography-top index on scalar-grid, used to vertically shift initial profiles
230    INTEGER(iwp) ::  nz_u_shift_l         !< topography-top index on u-grid, used to vertically shift initial profiles
231    INTEGER(iwp) ::  nz_v_shift_l         !< topography-top index on v-grid, used to vertically shift initial profiles
232    INTEGER(iwp) ::  nz_w_shift_l         !< topography-top index on w-grid, used to vertically shift initial profiles
233    INTEGER(iwp) ::  nz_s_shift_l         !< topography-top index on scalar-grid, used to vertically shift initial profiles
234    INTEGER(iwp) ::  nzt_l                !< index of top PE boundary for multigrid level
235    INTEGER(iwp) ::  pids_static_id       !< file id for static input file
236    INTEGER(iwp) ::  sr                   !< index of statistic region
237
238    INTEGER(iwp), DIMENSION(:), ALLOCATABLE   ::  ngp_2dh_l  !< toal number of horizontal grid points in statistical region on subdomain
239
240    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  ngp_2dh_outer_l    !< number of horizontal non-wall bounded grid points on subdomain
241    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  ngp_2dh_s_inner_l  !< number of horizontal non-topography grid points on subdomain
242
243
244   
245    REAL(wp), DIMENSION(:), ALLOCATABLE ::  init_l        !< dummy array used for averaging 3D data to obtain inital profiles
246    REAL(wp), DIMENSION(:), ALLOCATABLE ::  p_hydrostatic !< hydrostatic pressure
247
248    REAL(wp) ::  dx_l !< grid spacing along x on different multigrid level
249    REAL(wp) ::  dy_l !< grid spacing along y on different multigrid level
250
251    REAL(wp), DIMENSION(1:3) ::  volume_flow_area_l     !< area of lateral and top model domain surface on local subdomain
252    REAL(wp), DIMENSION(1:3) ::  volume_flow_initial_l  !< initial volume flow into model domain
253
254    REAL(wp), DIMENSION(:), ALLOCATABLE ::  mean_surface_level_height_l !< mean surface level height on subdomain
255    REAL(wp), DIMENSION(:), ALLOCATABLE ::  ngp_3d_inner_l    !< total number of non-topography grid points on subdomain
256    REAL(wp), DIMENSION(:), ALLOCATABLE ::  ngp_3d_inner_tmp  !< total number of non-topography grid points
257
258    TYPE(real_2d) ::  tmp_2d !< temporary variable to input additional surface-data from static file
259   
260    CALL location_message( 'model initialization', 'start' )
261
262    IF ( debug_output )  CALL debug_message( 'allocating arrays', 'start' )
263!
264!-- Allocate arrays
265    ALLOCATE( mean_surface_level_height(0:statistic_regions),                  &
266              mean_surface_level_height_l(0:statistic_regions),                &
267              ngp_2dh(0:statistic_regions), ngp_2dh_l(0:statistic_regions),    &
268              ngp_3d(0:statistic_regions),                                     &
269              ngp_3d_inner(0:statistic_regions),                               &
270              ngp_3d_inner_l(0:statistic_regions),                             &
271              ngp_3d_inner_tmp(0:statistic_regions),                           &
272              sums_divnew_l(0:statistic_regions),                              &
273              sums_divold_l(0:statistic_regions) )
274    ALLOCATE( dp_smooth_factor(nzb:nzt), rdf(nzb+1:nzt), rdf_sc(nzb+1:nzt) )
275    ALLOCATE( ngp_2dh_outer(nzb:nzt+1,0:statistic_regions),                    &
276              ngp_2dh_outer_l(nzb:nzt+1,0:statistic_regions),                  &
277              ngp_2dh_s_inner(nzb:nzt+1,0:statistic_regions),                  &
278              ngp_2dh_s_inner_l(nzb:nzt+1,0:statistic_regions),                &
279              rmask(nysg:nyng,nxlg:nxrg,0:statistic_regions),                  &
280              sums(nzb:nzt+1,pr_palm+max_pr_user+max_pr_cs+max_pr_salsa),      &
281              sums_l(nzb:nzt+1,pr_palm+max_pr_user+max_pr_cs+max_pr_salsa,0:threads_per_task-1),   &
282              sums_l_l(nzb:nzt+1,0:statistic_regions,0:threads_per_task-1),    &
283              sums_wsts_bc_l(nzb:nzt+1,0:statistic_regions) )
284    ALLOCATE( ts_value(dots_max,0:statistic_regions) )
285    ALLOCATE( ptdf_x(nxlg:nxrg), ptdf_y(nysg:nyng) )
286
287    ALLOCATE( d(nzb+1:nzt,nys:nyn,nxl:nxr),                                    &
288              p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                &
289              tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
290
291    ALLOCATE( pt_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
292              pt_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
293              u_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
294              u_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
295              u_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
296              v_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
297              v_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
298              v_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
299              w_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
300              w_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
301              w_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
302    IF (  .NOT.  neutral )  THEN
303       ALLOCATE( pt_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
304    ENDIF
305!
306!-- Pre-set masks for regional statistics. Default is the total model domain.
307!-- Ghost points are excluded because counting values at the ghost boundaries
308!-- would bias the statistics
309    rmask = 1.0_wp
310    rmask(:,nxlg:nxl-1,:) = 0.0_wp;  rmask(:,nxr+1:nxrg,:) = 0.0_wp
311    rmask(nysg:nys-1,:,:) = 0.0_wp;  rmask(nyn+1:nyng,:,:) = 0.0_wp
312!
313!-- Following array is required for perturbation pressure within the iterative
314!-- pressure solvers. For the multistep schemes (Runge-Kutta), array p holds
315!-- the weighted average of the substeps and cannot be used in the Poisson
316!-- solver.
317    IF ( psolver == 'sor' )  THEN
318       ALLOCATE( p_loc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
319    ELSEIF ( psolver(1:9) == 'multigrid' )  THEN
320!
321!--    For performance reasons, multigrid is using one ghost layer only
322       ALLOCATE( p_loc(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
323    ENDIF
324
325!
326!-- Array for storing constant coeffficients of the tridiagonal solver
327    IF ( psolver == 'poisfft' )  THEN
328       ALLOCATE( tri(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1,2) )
329       ALLOCATE( tric(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1) )
330    ENDIF
331
332    IF ( humidity )  THEN
333!
334!--    3D-humidity
335       ALLOCATE( q_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
336                 q_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
337                 q_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
338                 vpt_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 
339    ENDIF   
340   
341    IF ( passive_scalar )  THEN
342
343!
344!--    3D scalar arrays
345       ALLOCATE( s_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
346                 s_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
347                 s_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
348
349    ENDIF
350
351!
352!-- Allocate and set 1d-profiles for Stokes drift velocity. It may be set to
353!-- non-zero values later in ocean_init
354    ALLOCATE( u_stokes_zu(nzb:nzt+1), u_stokes_zw(nzb:nzt+1),                  &
355              v_stokes_zu(nzb:nzt+1), v_stokes_zw(nzb:nzt+1) )
356    u_stokes_zu(:) = 0.0_wp
357    u_stokes_zw(:) = 0.0_wp
358    v_stokes_zu(:) = 0.0_wp
359    v_stokes_zw(:) = 0.0_wp
360
361!
362!-- Allocation of anelastic and Boussinesq approximation specific arrays
363    ALLOCATE( p_hydrostatic(nzb:nzt+1) )
364    ALLOCATE( rho_air(nzb:nzt+1) )
365    ALLOCATE( rho_air_zw(nzb:nzt+1) )
366    ALLOCATE( drho_air(nzb:nzt+1) )
367    ALLOCATE( drho_air_zw(nzb:nzt+1) )
368!
369!-- Density profile calculation for anelastic and Boussinesq approximation
370!-- In case of a Boussinesq approximation, a constant density is calculated
371!-- mainly for output purposes. This density do not need to be considered
372!-- in the model's system of equations.
373    IF ( TRIM( approximation ) == 'anelastic' )  THEN
374       DO  k = nzb, nzt+1
375          p_hydrostatic(k) = barometric_formula(zu(k), pt_surface *            & 
376                             exner_function(surface_pressure * 100.0_wp),      &
377                             surface_pressure * 100.0_wp)
378         
379          rho_air(k) = ideal_gas_law_rho_pt(p_hydrostatic(k), pt_init(k))
380       ENDDO
381       
382       DO  k = nzb, nzt
383          rho_air_zw(k) = 0.5_wp * ( rho_air(k) + rho_air(k+1) )
384       ENDDO
385       
386       rho_air_zw(nzt+1)  = rho_air_zw(nzt)                                    &
387                            + 2.0_wp * ( rho_air(nzt+1) - rho_air_zw(nzt)  )
388                           
389    ELSE
390       DO  k = nzb, nzt+1
391          p_hydrostatic(k) = barometric_formula(zu(nzb), pt_surface *          & 
392                             exner_function(surface_pressure * 100.0_wp),      &
393                             surface_pressure * 100.0_wp)
394
395          rho_air(k) = ideal_gas_law_rho_pt(p_hydrostatic(k), pt_init(nzb))
396       ENDDO
397       
398       DO  k = nzb, nzt
399          rho_air_zw(k) = 0.5_wp * ( rho_air(k) + rho_air(k+1) )
400       ENDDO
401       
402       rho_air_zw(nzt+1)  = rho_air_zw(nzt)                                    &
403                            + 2.0_wp * ( rho_air(nzt+1) - rho_air_zw(nzt)  )
404                           
405    ENDIF
406!
407!-- compute the inverse density array in order to avoid expencive divisions
408    drho_air    = 1.0_wp / rho_air
409    drho_air_zw = 1.0_wp / rho_air_zw
410
411!
412!-- Allocation of flux conversion arrays
413    ALLOCATE( heatflux_input_conversion(nzb:nzt+1) )
414    ALLOCATE( waterflux_input_conversion(nzb:nzt+1) )
415    ALLOCATE( momentumflux_input_conversion(nzb:nzt+1) )
416    ALLOCATE( heatflux_output_conversion(nzb:nzt+1) )
417    ALLOCATE( waterflux_output_conversion(nzb:nzt+1) )
418    ALLOCATE( momentumflux_output_conversion(nzb:nzt+1) )
419
420!
421!-- calculate flux conversion factors according to approximation and in-/output mode
422    DO  k = nzb, nzt+1
423
424        IF ( TRIM( flux_input_mode ) == 'kinematic' )  THEN
425            heatflux_input_conversion(k)      = rho_air_zw(k)
426            waterflux_input_conversion(k)     = rho_air_zw(k)
427            momentumflux_input_conversion(k)  = rho_air_zw(k)
428        ELSEIF ( TRIM( flux_input_mode ) == 'dynamic' ) THEN
429            heatflux_input_conversion(k)      = 1.0_wp / c_p
430            waterflux_input_conversion(k)     = 1.0_wp / l_v
431            momentumflux_input_conversion(k)  = 1.0_wp
432        ENDIF
433
434        IF ( TRIM( flux_output_mode ) == 'kinematic' )  THEN
435            heatflux_output_conversion(k)     = drho_air_zw(k)
436            waterflux_output_conversion(k)    = drho_air_zw(k)
437            momentumflux_output_conversion(k) = drho_air_zw(k)
438        ELSEIF ( TRIM( flux_output_mode ) == 'dynamic' ) THEN
439            heatflux_output_conversion(k)     = c_p
440            waterflux_output_conversion(k)    = l_v
441            momentumflux_output_conversion(k) = 1.0_wp
442        ENDIF
443
444        IF ( .NOT. humidity ) THEN
445            waterflux_input_conversion(k)  = 1.0_wp
446            waterflux_output_conversion(k) = 1.0_wp
447        ENDIF
448
449    ENDDO
450
451!
452!-- In case of multigrid method, compute grid lengths and grid factors for the
453!-- grid levels with respective density on each grid
454    IF ( psolver(1:9) == 'multigrid' )  THEN
455
456       ALLOCATE( ddx2_mg(maximum_grid_level) )
457       ALLOCATE( ddy2_mg(maximum_grid_level) )
458       ALLOCATE( dzu_mg(nzb+1:nzt+1,maximum_grid_level) )
459       ALLOCATE( dzw_mg(nzb+1:nzt+1,maximum_grid_level) )
460       ALLOCATE( f1_mg(nzb+1:nzt,maximum_grid_level) )
461       ALLOCATE( f2_mg(nzb+1:nzt,maximum_grid_level) )
462       ALLOCATE( f3_mg(nzb+1:nzt,maximum_grid_level) )
463       ALLOCATE( rho_air_mg(nzb:nzt+1,maximum_grid_level) )
464       ALLOCATE( rho_air_zw_mg(nzb:nzt+1,maximum_grid_level) )
465
466       dzu_mg(:,maximum_grid_level) = dzu
467       rho_air_mg(:,maximum_grid_level) = rho_air
468!       
469!--    Next line to ensure an equally spaced grid.
470       dzu_mg(1,maximum_grid_level) = dzu(2)
471       rho_air_mg(nzb,maximum_grid_level) = rho_air(nzb) +                     &
472                                             (rho_air(nzb) - rho_air(nzb+1))
473
474       dzw_mg(:,maximum_grid_level) = dzw
475       rho_air_zw_mg(:,maximum_grid_level) = rho_air_zw
476       nzt_l = nzt
477       DO  l = maximum_grid_level-1, 1, -1
478           dzu_mg(nzb+1,l) = 2.0_wp * dzu_mg(nzb+1,l+1)
479           dzw_mg(nzb+1,l) = 2.0_wp * dzw_mg(nzb+1,l+1)
480           rho_air_mg(nzb,l)    = rho_air_mg(nzb,l+1) + (rho_air_mg(nzb,l+1) - rho_air_mg(nzb+1,l+1))
481           rho_air_zw_mg(nzb,l) = rho_air_zw_mg(nzb,l+1) + (rho_air_zw_mg(nzb,l+1) - rho_air_zw_mg(nzb+1,l+1))
482           rho_air_mg(nzb+1,l)    = rho_air_mg(nzb+1,l+1)
483           rho_air_zw_mg(nzb+1,l) = rho_air_zw_mg(nzb+1,l+1)
484           nzt_l = nzt_l / 2
485           DO  k = 2, nzt_l+1
486              dzu_mg(k,l) = dzu_mg(2*k-2,l+1) + dzu_mg(2*k-1,l+1)
487              dzw_mg(k,l) = dzw_mg(2*k-2,l+1) + dzw_mg(2*k-1,l+1)
488              rho_air_mg(k,l)    = rho_air_mg(2*k-1,l+1)
489              rho_air_zw_mg(k,l) = rho_air_zw_mg(2*k-1,l+1)
490           ENDDO
491       ENDDO
492
493       nzt_l = nzt
494       dx_l  = dx
495       dy_l  = dy
496       DO  l = maximum_grid_level, 1, -1
497          ddx2_mg(l) = 1.0_wp / dx_l**2
498          ddy2_mg(l) = 1.0_wp / dy_l**2
499          DO  k = nzb+1, nzt_l
500             f2_mg(k,l) = rho_air_zw_mg(k,l) / ( dzu_mg(k+1,l) * dzw_mg(k,l) )
501             f3_mg(k,l) = rho_air_zw_mg(k-1,l) / ( dzu_mg(k,l)   * dzw_mg(k,l) )
502             f1_mg(k,l) = 2.0_wp * ( ddx2_mg(l) + ddy2_mg(l) ) &
503                          * rho_air_mg(k,l) + f2_mg(k,l) + f3_mg(k,l)
504          ENDDO
505          nzt_l = nzt_l / 2
506          dx_l  = dx_l * 2.0_wp
507          dy_l  = dy_l * 2.0_wp
508       ENDDO
509
510    ENDIF
511
512!
513!-- 1D-array for large scale subsidence velocity
514    IF ( .NOT. ALLOCATED( w_subs ) )  THEN
515       ALLOCATE ( w_subs(nzb:nzt+1) )
516       w_subs = 0.0_wp
517    ENDIF
518
519!
520!-- Arrays to store velocity data from t-dt and the phase speeds which
521!-- are needed for radiation boundary conditions
522    IF ( bc_radiation_l )  THEN
523       ALLOCATE( u_m_l(nzb:nzt+1,nysg:nyng,1:2),                               &
524                 v_m_l(nzb:nzt+1,nysg:nyng,0:1),                               &
525                 w_m_l(nzb:nzt+1,nysg:nyng,0:1) )
526    ENDIF
527    IF ( bc_radiation_r )  THEN
528       ALLOCATE( u_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx),                           &
529                 v_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx),                           &
530                 w_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx) )
531    ENDIF
532    IF ( bc_radiation_l  .OR.  bc_radiation_r )  THEN
533       ALLOCATE( c_u(nzb:nzt+1,nysg:nyng), c_v(nzb:nzt+1,nysg:nyng),           &
534                 c_w(nzb:nzt+1,nysg:nyng) )
535    ENDIF
536    IF ( bc_radiation_s )  THEN
537       ALLOCATE( u_m_s(nzb:nzt+1,0:1,nxlg:nxrg),                               &
538                 v_m_s(nzb:nzt+1,1:2,nxlg:nxrg),                               &
539                 w_m_s(nzb:nzt+1,0:1,nxlg:nxrg) )
540    ENDIF
541    IF ( bc_radiation_n )  THEN
542       ALLOCATE( u_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg),                           &
543                 v_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg),                           &
544                 w_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg) )
545    ENDIF
546    IF ( bc_radiation_s  .OR.  bc_radiation_n )  THEN
547       ALLOCATE( c_u(nzb:nzt+1,nxlg:nxrg), c_v(nzb:nzt+1,nxlg:nxrg),           &
548                 c_w(nzb:nzt+1,nxlg:nxrg) )
549    ENDIF
550    IF ( bc_radiation_l  .OR.  bc_radiation_r  .OR.  bc_radiation_s  .OR.      &
551         bc_radiation_n )  THEN
552       ALLOCATE( c_u_m_l(nzb:nzt+1), c_v_m_l(nzb:nzt+1), c_w_m_l(nzb:nzt+1) )                   
553       ALLOCATE( c_u_m(nzb:nzt+1), c_v_m(nzb:nzt+1), c_w_m(nzb:nzt+1) )
554    ENDIF
555
556!
557!-- Initial assignment of the pointers
558    IF ( .NOT. neutral )  THEN
559       pt => pt_1;  pt_p => pt_2;  tpt_m => pt_3
560    ELSE
561       pt => pt_1;  pt_p => pt_1;  tpt_m => pt_3
562    ENDIF
563    u  => u_1;   u_p  => u_2;   tu_m  => u_3
564    v  => v_1;   v_p  => v_2;   tv_m  => v_3
565    w  => w_1;   w_p  => w_2;   tw_m  => w_3
566
567    IF ( humidity )  THEN
568       q => q_1;  q_p => q_2;  tq_m => q_3
569       vpt  => vpt_1
570    ENDIF
571   
572    IF ( passive_scalar )  THEN
573       s => s_1;  s_p => s_2;  ts_m => s_3
574    ENDIF   
575
576!
577!-- Initialize surface arrays
578    CALL init_surface_arrays
579!
580!-- Allocate arrays for other modules
581    CALL module_interface_init_arrays
582
583
584!
585!-- Allocate arrays containing the RK coefficient for calculation of
586!-- perturbation pressure and turbulent fluxes. At this point values are
587!-- set for pressure calculation during initialization (where no timestep
588!-- is done). Further below the values needed within the timestep scheme
589!-- will be set.
590    ALLOCATE( weight_substep(1:intermediate_timestep_count_max),               &
591              weight_pres(1:intermediate_timestep_count_max) )
592    weight_substep = 1.0_wp
593    weight_pres    = 1.0_wp
594    intermediate_timestep_count = 0  ! needed when simulated_time = 0.0
595       
596    IF ( debug_output )  CALL debug_message( 'allocating arrays', 'end' )
597
598!
599!-- Initialize time series
600    ts_value = 0.0_wp
601
602!
603!-- Initialize local summation arrays for routine flow_statistics.
604!-- This is necessary because they may not yet have been initialized when they
605!-- are called from flow_statistics (or - depending on the chosen model run -
606!-- are never initialized)
607    sums_divnew_l      = 0.0_wp
608    sums_divold_l      = 0.0_wp
609    sums_l_l           = 0.0_wp
610    sums_wsts_bc_l     = 0.0_wp
611   
612!
613!-- Initialize model variables
614    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
615         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
616!
617!--    Initialization with provided input data derived from larger-scale model
618       IF ( INDEX( initializing_actions, 'inifor' ) /= 0 )  THEN
619          IF ( debug_output )  CALL debug_message( 'initializing with INIFOR', 'start' )
620!
621!--       Read initial 1D profiles or 3D data from NetCDF file, depending
622!--       on the provided level-of-detail.
623!--       At the moment, only u, v, w, pt and q are provided.
624          CALL netcdf_data_input_init_3d
625!
626!--       Please note, Inifor provides data from nzb+1 to nzt.
627!--       Bottom and top boundary conditions for Inifor profiles are already
628!--       set (just after reading), so that this is not necessary here.
629!--       Depending on the provided level-of-detail, initial Inifor data is
630!--       either stored on data type (lod=1), or directly on 3D arrays (lod=2).
631!--       In order to obtain also initial profiles in case of lod=2 (which
632!--       is required for e.g. damping), average over 3D data.
633          IF( init_3d%lod_u == 1 )  THEN
634             u_init = init_3d%u_init
635          ELSEIF( init_3d%lod_u == 2 )  THEN
636             ALLOCATE( init_l(nzb:nzt+1) ) 
637             DO  k = nzb, nzt+1
638                init_l(k) = SUM( u(k,nys:nyn,nxl:nxr) )
639             ENDDO
640             init_l = init_l / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
641
642#if defined( __parallel )
643             CALL MPI_ALLREDUCE( init_l, u_init, nzt+1-nzb+1,                  &
644                                 MPI_REAL, MPI_SUM, comm2d, ierr )
645#else
646             u_init = init_l
647#endif
648             DEALLOCATE( init_l )
649
650          ENDIF
651           
652          IF( init_3d%lod_v == 1 )  THEN 
653             v_init = init_3d%v_init
654          ELSEIF( init_3d%lod_v == 2 )  THEN
655             ALLOCATE( init_l(nzb:nzt+1) ) 
656             DO  k = nzb, nzt+1
657                init_l(k) = SUM( v(k,nys:nyn,nxl:nxr) )
658             ENDDO
659             init_l = init_l / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
660
661#if defined( __parallel )
662             CALL MPI_ALLREDUCE( init_l, v_init, nzt+1-nzb+1,                  &
663                                 MPI_REAL, MPI_SUM, comm2d, ierr )
664#else
665             v_init = init_l
666#endif
667             DEALLOCATE( init_l )
668          ENDIF
669          IF( .NOT. neutral )  THEN
670             IF( init_3d%lod_pt == 1 )  THEN
671                pt_init = init_3d%pt_init
672             ELSEIF( init_3d%lod_pt == 2 )  THEN
673                ALLOCATE( init_l(nzb:nzt+1) ) 
674                DO  k = nzb, nzt+1
675                   init_l(k) = SUM( pt(k,nys:nyn,nxl:nxr) )
676                ENDDO
677                init_l = init_l / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
678
679#if defined( __parallel )
680                CALL MPI_ALLREDUCE( init_l, pt_init, nzt+1-nzb+1,               &
681                                    MPI_REAL, MPI_SUM, comm2d, ierr )
682#else
683                pt_init = init_l
684#endif
685                DEALLOCATE( init_l )
686             ENDIF
687          ENDIF
688
689
690          IF( humidity )  THEN
691             IF( init_3d%lod_q == 1 )  THEN
692                q_init = init_3d%q_init
693             ELSEIF( init_3d%lod_q == 2 )  THEN
694                ALLOCATE( init_l(nzb:nzt+1) ) 
695                DO  k = nzb, nzt+1
696                   init_l(k) = SUM( q(k,nys:nyn,nxl:nxr) )
697                ENDDO
698                init_l = init_l / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
699
700#if defined( __parallel )
701                CALL MPI_ALLREDUCE( init_l, q_init, nzt+1-nzb+1,               &
702                                    MPI_REAL, MPI_SUM, comm2d, ierr )
703#else
704                q_init = init_l
705#endif
706                DEALLOCATE( init_l )
707             ENDIF
708          ENDIF
709
710!
711!--       Write initial profiles onto 3D arrays.
712!--       Work-around, 3D initialization of u,v,w creates artificial
713!--       structures wich correlate with the processor grid. The reason
714!--       for this is still unknown. To work-around this, 3D initialization
715!--       will be effectively reduce to a 1D initialization where no such
716!--       artificial structures appear.
717          DO  i = nxlg, nxrg
718             DO  j = nysg, nyng
719                IF( init_3d%lod_u == 1  .OR.  init_3d%lod_u == 2 )             &
720                   u(:,j,i) = u_init(:)
721                IF( init_3d%lod_v == 1  .OR.  init_3d%lod_u == 2 )             &
722                   v(:,j,i) = v_init(:)
723                IF( .NOT. neutral  .AND.                                       &
724                    ( init_3d%lod_pt == 1  .OR.  init_3d%lod_pt == 2 ) )       &
725                   pt(:,j,i) = pt_init(:)
726                IF( humidity  .AND.                                            &
727                    ( init_3d%lod_q == 1  .OR.  init_3d%lod_q == 2 ) )         &
728                   q(:,j,i) = q_init(:)
729             ENDDO
730          ENDDO
731!
732!--       Set geostrophic wind components. 
733          IF ( init_3d%from_file_ug )  THEN
734             ug(:) = init_3d%ug_init(:)
735          ENDIF
736          IF ( init_3d%from_file_vg )  THEN
737             vg(:) = init_3d%vg_init(:)
738          ENDIF
739!
740!--       Set bottom and top boundary condition for geostrophic wind
741          ug(nzt+1) = ug(nzt)
742          vg(nzt+1) = vg(nzt)
743          ug(nzb)   = ug(nzb+1)
744          vg(nzb)   = vg(nzb+1)
745!
746!--       Set inital w to 0
747          w = 0.0_wp
748
749          IF ( passive_scalar )  THEN
750             DO  i = nxlg, nxrg
751                DO  j = nysg, nyng
752                   s(:,j,i) = s_init
753                ENDDO
754             ENDDO
755          ENDIF
756
757!
758!--       Set velocity components at non-atmospheric / oceanic grid points to
759!--       zero.
760          u = MERGE( u, 0.0_wp, BTEST( wall_flags_0, 1 ) )
761          v = MERGE( v, 0.0_wp, BTEST( wall_flags_0, 2 ) )
762          w = MERGE( w, 0.0_wp, BTEST( wall_flags_0, 3 ) )
763!
764!--       Initialize surface variables, e.g. friction velocity, momentum
765!--       fluxes, etc.
766          CALL init_surfaces
767
768          IF ( debug_output )  CALL debug_message( 'initializing with INIFOR', 'end' )
769!
770!--    Initialization via computed 1D-model profiles
771       ELSEIF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
772
773          IF ( debug_output )  CALL debug_message( 'initializing with 1D model profiles', 'start' )
774!
775!--       Use solutions of the 1D model as initial profiles,
776!--       start 1D model
777          CALL init_1d_model
778!
779!--       Transfer initial profiles to the arrays of the 3D model
780          DO  i = nxlg, nxrg
781             DO  j = nysg, nyng
782                pt(:,j,i) = pt_init
783                u(:,j,i)  = u1d
784                v(:,j,i)  = v1d
785             ENDDO
786          ENDDO
787
788          IF ( humidity )  THEN
789             DO  i = nxlg, nxrg
790                DO  j = nysg, nyng
791                   q(:,j,i) = q_init
792                ENDDO
793             ENDDO
794          ENDIF
795
796          IF ( passive_scalar )  THEN
797             DO  i = nxlg, nxrg
798                DO  j = nysg, nyng
799                   s(:,j,i) = s_init
800                ENDDO
801             ENDDO   
802          ENDIF
803!
804!--          Store initial profiles for output purposes etc.
805          IF ( .NOT. constant_diffusion )  THEN
806             hom(:,1,25,:) = SPREAD( l1d, 2, statistic_regions+1 )
807          ENDIF
808!
809!--       Set velocities back to zero
810          u = MERGE( u, 0.0_wp, BTEST( wall_flags_0, 1 ) )
811          v = MERGE( v, 0.0_wp, BTEST( wall_flags_0, 2 ) )         
812!
813!--       WARNING: The extra boundary conditions set after running the
814!--       -------  1D model impose an error on the divergence one layer
815!--                below the topography; need to correct later
816!--       ATTENTION: Provisional correction for Piacsek & Williams
817!--       ---------  advection scheme: keep u and v zero one layer below
818!--                  the topography.
819          IF ( ibc_uv_b == 1 )  THEN
820!
821!--          Neumann condition
822             DO  i = nxl-1, nxr+1
823                DO  j = nys-1, nyn+1
824                   u(nzb,j,i) = u(nzb+1,j,i)
825                   v(nzb,j,i) = v(nzb+1,j,i)
826                ENDDO
827             ENDDO
828
829          ENDIF
830!
831!--       Initialize surface variables, e.g. friction velocity, momentum
832!--       fluxes, etc.
833          CALL init_surfaces
834
835          IF ( debug_output )  CALL debug_message( 'initializing with 1D model profiles', 'end' )
836
837       ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 )    &
838       THEN
839
840          IF ( debug_output )  CALL debug_message( 'initializing with constant profiles', 'start' )
841
842!
843!--       Use constructed initial profiles (velocity constant with height,
844!--       temperature profile with constant gradient)
845          DO  i = nxlg, nxrg
846             DO  j = nysg, nyng
847                pt(:,j,i) = pt_init
848                u(:,j,i)  = u_init
849                v(:,j,i)  = v_init
850             ENDDO
851          ENDDO
852!
853!--       Mask topography
854          u = MERGE( u, 0.0_wp, BTEST( wall_flags_0, 1 ) )
855          v = MERGE( v, 0.0_wp, BTEST( wall_flags_0, 2 ) )
856!
857!--       Set initial horizontal velocities at the lowest computational grid
858!--       levels to zero in order to avoid too small time steps caused by the
859!--       diffusion limit in the initial phase of a run (at k=1, dz/2 occurs
860!--       in the limiting formula!).
861!--       Please note, in case land- or urban-surface model is used and a
862!--       spinup is applied, masking the lowest computational level is not
863!--       possible as MOST as well as energy-balance parametrizations will not
864!--       work with zero wind velocity.
865          IF ( ibc_uv_b /= 1  .AND.  .NOT.  spinup )  THEN
866             DO  i = nxlg, nxrg
867                DO  j = nysg, nyng
868                   DO  k = nzb, nzt
869                      u(k,j,i) = MERGE( u(k,j,i), 0.0_wp,                      &
870                                        BTEST( wall_flags_0(k,j,i), 20 ) )
871                      v(k,j,i) = MERGE( v(k,j,i), 0.0_wp,                      &
872                                        BTEST( wall_flags_0(k,j,i), 21 ) )
873                   ENDDO
874                ENDDO
875             ENDDO
876          ENDIF
877
878          IF ( humidity )  THEN
879             DO  i = nxlg, nxrg
880                DO  j = nysg, nyng
881                   q(:,j,i) = q_init
882                ENDDO
883             ENDDO
884          ENDIF
885         
886          IF ( passive_scalar )  THEN
887             DO  i = nxlg, nxrg
888                DO  j = nysg, nyng
889                   s(:,j,i) = s_init
890                ENDDO
891             ENDDO
892          ENDIF
893
894!
895!--       Compute initial temperature field and other constants used in case
896!--       of a sloping surface
897          IF ( sloping_surface )  CALL init_slope
898!
899!--       Initialize surface variables, e.g. friction velocity, momentum
900!--       fluxes, etc.
901          CALL init_surfaces
902         
903          IF ( debug_output )  CALL debug_message( 'initializing with constant profiles', 'end' )
904
905       ELSEIF ( INDEX(initializing_actions, 'by_user') /= 0 )                  &
906       THEN
907
908          IF ( debug_output )  CALL debug_message( 'initializing by user', 'start' )
909!
910!--       Pre-initialize surface variables, i.e. setting start- and end-indices
911!--       at each (j,i)-location. Please note, this does not supersede
912!--       user-defined initialization of surface quantities.
913          CALL init_surfaces
914!
915!--       Initialization will completely be done by the user
916          CALL user_init_3d_model
917
918          IF ( debug_output )  CALL debug_message( 'initializing by user', 'end' )
919
920       ENDIF
921
922       IF ( debug_output )  CALL debug_message( 'initializing statistics, boundary conditions, etc.', 'start' )
923
924!
925!--    Bottom boundary
926       IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2  )  THEN
927          u(nzb,:,:) = 0.0_wp
928          v(nzb,:,:) = 0.0_wp
929       ENDIF
930
931!
932!--    Apply channel flow boundary condition
933       IF ( TRIM( bc_uv_t ) == 'dirichlet_0' )  THEN
934          u(nzt+1,:,:) = 0.0_wp
935          v(nzt+1,:,:) = 0.0_wp
936       ENDIF
937
938!
939!--    Calculate virtual potential temperature
940       IF ( humidity )  vpt = pt * ( 1.0_wp + 0.61_wp * q )
941
942!
943!--    Store initial profiles for output purposes etc.. Please note, in case of
944!--    initialization of u, v, w, pt, and q via output data derived from larger
945!--    scale models, data will not be horizontally homogeneous. Actually, a mean
946!--    profile should be calculated before.   
947       hom(:,1,5,:) = SPREAD( u(:,nys,nxl), 2, statistic_regions+1 )
948       hom(:,1,6,:) = SPREAD( v(:,nys,nxl), 2, statistic_regions+1 )
949       IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2)  THEN
950          hom(nzb,1,5,:) = 0.0_wp
951          hom(nzb,1,6,:) = 0.0_wp
952       ENDIF
953       hom(:,1,7,:)  = SPREAD( pt(:,nys,nxl), 2, statistic_regions+1 )
954
955       IF ( humidity )  THEN
956!
957!--       Store initial profile of total water content, virtual potential
958!--       temperature
959          hom(:,1,26,:) = SPREAD(   q(:,nys,nxl), 2, statistic_regions+1 )
960          hom(:,1,29,:) = SPREAD( vpt(:,nys,nxl), 2, statistic_regions+1 )
961!
962!--       Store initial profile of mixing ratio and potential
963!--       temperature
964          IF ( bulk_cloud_model  .OR.  cloud_droplets ) THEN
965             hom(:,1,27,:) = SPREAD(  q(:,nys,nxl), 2, statistic_regions+1 )
966             hom(:,1,28,:) = SPREAD( pt(:,nys,nxl), 2, statistic_regions+1 )
967          ENDIF
968       ENDIF
969
970!
971!--    Store initial scalar profile
972       IF ( passive_scalar )  THEN
973          hom(:,1,121,:) = SPREAD(  s(:,nys,nxl), 2, statistic_regions+1 )
974       ENDIF
975
976!
977!--    Initialize the random number generators (from numerical recipes)
978       CALL random_function_ini
979       
980       IF ( random_generator == 'random-parallel' )  THEN
981          CALL init_parallel_random_generator( nx, nys, nyn, nxl, nxr )
982       ENDIF
983!
984!--    Set the reference state to be used in the buoyancy terms (for ocean runs
985!--    the reference state will be set (overwritten) in init_ocean)
986       IF ( use_single_reference_value )  THEN
987          IF (  .NOT.  humidity )  THEN
988             ref_state(:) = pt_reference
989          ELSE
990             ref_state(:) = vpt_reference
991          ENDIF
992       ELSE
993          IF (  .NOT.  humidity )  THEN
994             ref_state(:) = pt_init(:)
995          ELSE
996             ref_state(:) = vpt(:,nys,nxl)
997          ENDIF
998       ENDIF
999
1000!
1001!--    For the moment, vertical velocity is zero
1002       w = 0.0_wp
1003
1004!
1005!--    Initialize array sums (must be defined in first call of pres)
1006       sums = 0.0_wp
1007
1008!
1009!--    In case of iterative solvers, p must get an initial value
1010       IF ( psolver(1:9) == 'multigrid'  .OR.  psolver == 'sor' )  p = 0.0_wp
1011!
1012!--    Impose vortex with vertical axis on the initial velocity profile
1013       IF ( INDEX( initializing_actions, 'initialize_vortex' ) /= 0 )  THEN
1014          CALL init_rankine
1015       ENDIF
1016
1017!
1018!--    Impose temperature anomaly (advection test only) or warm air bubble
1019!--    close to surface
1020       IF ( INDEX( initializing_actions, 'initialize_ptanom' ) /= 0  .OR.  &
1021            INDEX( initializing_actions, 'initialize_bubble' ) /= 0  )  THEN
1022          CALL init_pt_anomaly
1023       ENDIF
1024       
1025!
1026!--    If required, change the surface temperature at the start of the 3D run
1027       IF ( pt_surface_initial_change /= 0.0_wp )  THEN
1028          pt(nzb,:,:) = pt(nzb,:,:) + pt_surface_initial_change
1029       ENDIF
1030
1031!
1032!--    If required, change the surface humidity/scalar at the start of the 3D
1033!--    run
1034       IF ( humidity  .AND.  q_surface_initial_change /= 0.0_wp )              &
1035          q(nzb,:,:) = q(nzb,:,:) + q_surface_initial_change
1036         
1037       IF ( passive_scalar .AND.  s_surface_initial_change /= 0.0_wp )         &
1038          s(nzb,:,:) = s(nzb,:,:) + s_surface_initial_change
1039       
1040
1041!
1042!--    Initialize old and new time levels.
1043       tpt_m = 0.0_wp; tu_m = 0.0_wp; tv_m = 0.0_wp; tw_m = 0.0_wp
1044       pt_p = pt; u_p = u; v_p = v; w_p = w
1045
1046       IF ( humidity  )  THEN
1047          tq_m = 0.0_wp
1048          q_p = q
1049       ENDIF
1050       
1051       IF ( passive_scalar )  THEN
1052          ts_m = 0.0_wp
1053          s_p  = s
1054       ENDIF       
1055
1056       IF ( debug_output )  CALL debug_message( 'initializing statistics, boundary conditions, etc.', 'end' )
1057
1058    ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data'  .OR.         &
1059             TRIM( initializing_actions ) == 'cyclic_fill' )                   &
1060    THEN
1061
1062       IF ( debug_output )  CALL debug_message( 'initializing in case of restart / cyclic_fill', 'start' )
1063!
1064!--    Initialize surface elements and its attributes, e.g. heat- and
1065!--    momentumfluxes, roughness, scaling parameters. As number of surface
1066!--    elements might be different between runs, e.g. in case of cyclic fill,
1067!--    and not all surface elements are read, surface elements need to be
1068!--    initialized before.
1069!--    Please note, in case of cyclic fill, surfaces should be initialized
1070!--    after restart data is read, else, individual settings of surface
1071!--    parameters will be overwritten from data of precursor run, hence,
1072!--    init_surfaces is called a second time after reading the restart data.
1073       CALL init_surfaces                       
1074!
1075!--    When reading data for cyclic fill of 3D prerun data files, read
1076!--    some of the global variables from the restart file which are required
1077!--    for initializing the inflow
1078       IF ( TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
1079
1080          DO  i = 0, io_blocks-1
1081             IF ( i == io_group )  THEN
1082                CALL rrd_read_parts_of_global
1083             ENDIF
1084#if defined( __parallel )
1085             CALL MPI_BARRIER( comm2d, ierr )
1086#endif
1087          ENDDO
1088
1089       ENDIF
1090
1091!
1092!--    Read processor specific binary data from restart file
1093       DO  i = 0, io_blocks-1
1094          IF ( i == io_group )  THEN
1095             CALL rrd_local
1096          ENDIF
1097#if defined( __parallel )
1098          CALL MPI_BARRIER( comm2d, ierr )
1099#endif
1100       ENDDO
1101!
1102!--    In case of cyclic fill, call init_surfaces a second time, so that
1103!--    surface properties such as heat fluxes are initialized as prescribed.
1104       IF ( TRIM( initializing_actions ) == 'cyclic_fill' )                    &
1105          CALL init_surfaces
1106
1107!
1108!--    In case of complex terrain and cyclic fill method as initialization,
1109!--    shift initial data in the vertical direction for each point in the
1110!--    x-y-plane depending on local surface height
1111       IF ( complex_terrain  .AND.                                             &
1112            TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
1113          DO  i = nxlg, nxrg
1114             DO  j = nysg, nyng
1115                nz_u_shift = topo_top_ind(j,i,1)
1116                nz_v_shift = topo_top_ind(j,i,2)
1117                nz_w_shift = topo_top_ind(j,i,3)
1118                nz_s_shift = topo_top_ind(j,i,0)
1119
1120                u(nz_u_shift:nzt+1,j,i)  = u(0:nzt+1-nz_u_shift,j,i)               
1121
1122                v(nz_v_shift:nzt+1,j,i)  = v(0:nzt+1-nz_v_shift,j,i)
1123
1124                w(nz_w_shift:nzt+1,j,i)  = w(0:nzt+1-nz_w_shift,j,i)
1125
1126                p(nz_s_shift:nzt+1,j,i)  =  p(0:nzt+1-nz_s_shift,j,i)
1127                pt(nz_s_shift:nzt+1,j,i) = pt(0:nzt+1-nz_s_shift,j,i)
1128             ENDDO
1129          ENDDO
1130       ENDIF
1131!
1132!--    Initialization of the turbulence recycling method
1133       IF ( TRIM( initializing_actions ) == 'cyclic_fill'  .AND.               &
1134            turbulent_inflow )  THEN
1135!
1136!--       First store the profiles to be used at the inflow.
1137!--       These profiles are the (temporally) and horizontally averaged vertical
1138!--       profiles from the prerun. Alternatively, prescribed profiles
1139!--       for u,v-components can be used.
1140!--       Overwrite u_init, v_init, pt_init, q_init and s_init with the
1141!--       (temporally) and horizontally averaged vertical profiles from the end
1142!--       of the prerun, because these profiles shall be used as the basic state
1143!--       for the rayleigh damping and the pt_damping.
1144          ALLOCATE( mean_inflow_profiles(nzb:nzt+1,1:num_mean_inflow_profiles) )
1145
1146          IF ( use_prescribed_profile_data )  THEN
1147             mean_inflow_profiles(:,1) = u_init            ! u
1148             mean_inflow_profiles(:,2) = v_init            ! v
1149          ELSE
1150             mean_inflow_profiles(:,1) = hom_sum(:,1,0)    ! u
1151             u_init(:)                 = hom_sum(:,1,0)
1152             mean_inflow_profiles(:,2) = hom_sum(:,2,0)    ! v
1153             v_init(:)                 = hom_sum(:,2,0)
1154          ENDIF
1155          mean_inflow_profiles(:,4) = hom_sum(:,4,0)       ! pt
1156          pt_init(:)                = hom_sum(:,4,0)
1157          IF ( humidity )                                                      &
1158             mean_inflow_profiles(:,6) = hom_sum(:,41,0)   ! q
1159             q_init(:)                 = hom_sum(:,41,0)   
1160          IF ( passive_scalar )                                                &
1161             mean_inflow_profiles(:,7) = hom_sum(:,115,0)   ! s
1162             s_init(:)                 = hom_sum(:,115,0)
1163!
1164!--       In case of complex terrain, determine vertical displacement at inflow
1165!--       boundary and adjust mean inflow profiles
1166          IF ( complex_terrain )  THEN
1167             IF ( nxlg <= 0 .AND. nxrg >= 0 .AND. nysg <= 0 .AND. nyng >= 0 )  THEN
1168                nz_u_shift_l = topo_top_ind(j,i,1)
1169                nz_v_shift_l = topo_top_ind(j,i,2)
1170                nz_w_shift_l = topo_top_ind(j,i,3)
1171                nz_s_shift_l = topo_top_ind(j,i,0)
1172             ELSE
1173                nz_u_shift_l = 0
1174                nz_v_shift_l = 0
1175                nz_w_shift_l = 0
1176                nz_s_shift_l = 0
1177             ENDIF
1178
1179#if defined( __parallel )
1180             CALL MPI_ALLREDUCE(nz_u_shift_l, nz_u_shift, 1, MPI_INTEGER,      &
1181                                MPI_MAX, comm2d, ierr)
1182             CALL MPI_ALLREDUCE(nz_v_shift_l, nz_v_shift, 1, MPI_INTEGER,      &
1183                                MPI_MAX, comm2d, ierr)
1184             CALL MPI_ALLREDUCE(nz_w_shift_l, nz_w_shift, 1, MPI_INTEGER,      & 
1185                                MPI_MAX, comm2d, ierr)
1186             CALL MPI_ALLREDUCE(nz_s_shift_l, nz_s_shift, 1, MPI_INTEGER,      &
1187                                MPI_MAX, comm2d, ierr)
1188#else
1189             nz_u_shift = nz_u_shift_l
1190             nz_v_shift = nz_v_shift_l
1191             nz_w_shift = nz_w_shift_l
1192             nz_s_shift = nz_s_shift_l
1193#endif
1194
1195             mean_inflow_profiles(:,1) = 0.0_wp
1196             mean_inflow_profiles(nz_u_shift:nzt+1,1) = hom_sum(0:nzt+1-nz_u_shift,1,0)  ! u
1197
1198             mean_inflow_profiles(:,2) = 0.0_wp
1199             mean_inflow_profiles(nz_v_shift:nzt+1,2) = hom_sum(0:nzt+1-nz_v_shift,2,0)  ! v
1200
1201             mean_inflow_profiles(nz_s_shift:nzt+1,4) = hom_sum(0:nzt+1-nz_s_shift,4,0)  ! pt
1202
1203          ENDIF
1204
1205!
1206!--       If necessary, adjust the horizontal flow field to the prescribed
1207!--       profiles
1208          IF ( use_prescribed_profile_data )  THEN
1209             DO  i = nxlg, nxrg
1210                DO  j = nysg, nyng
1211                   DO  k = nzb, nzt+1
1212                      u(k,j,i) = u(k,j,i) - hom_sum(k,1,0) + u_init(k)
1213                      v(k,j,i) = v(k,j,i) - hom_sum(k,2,0) + v_init(k)
1214                   ENDDO
1215                ENDDO
1216             ENDDO
1217          ENDIF
1218
1219!
1220!--       Use these mean profiles at the inflow (provided that Dirichlet
1221!--       conditions are used)
1222          IF ( bc_dirichlet_l )  THEN
1223             DO  j = nysg, nyng
1224                DO  k = nzb, nzt+1
1225                   u(k,j,nxlg:-1)  = mean_inflow_profiles(k,1)
1226                   v(k,j,nxlg:-1)  = mean_inflow_profiles(k,2)
1227                   w(k,j,nxlg:-1)  = 0.0_wp
1228                   pt(k,j,nxlg:-1) = mean_inflow_profiles(k,4)
1229                   IF ( humidity )                                             &
1230                      q(k,j,nxlg:-1)  = mean_inflow_profiles(k,6)
1231                   IF ( passive_scalar )                                       &
1232                      s(k,j,nxlg:-1)  = mean_inflow_profiles(k,7)                     
1233                ENDDO
1234             ENDDO
1235          ENDIF
1236
1237!
1238!--       Calculate the damping factors to be used at the inflow. For a
1239!--       turbulent inflow the turbulent fluctuations have to be limited
1240!--       vertically because otherwise the turbulent inflow layer will grow
1241!--       in time.
1242          IF ( inflow_damping_height == 9999999.9_wp )  THEN
1243!
1244!--          Default: use the inversion height calculated by the prerun; if
1245!--          this is zero, inflow_damping_height must be explicitly
1246!--          specified.
1247             IF ( hom_sum(nzb+6,pr_palm,0) /= 0.0_wp )  THEN
1248                inflow_damping_height = hom_sum(nzb+6,pr_palm,0)
1249             ELSE
1250                WRITE( message_string, * ) 'inflow_damping_height must be ',   &
1251                     'explicitly specified because&the inversion height ',     &
1252                     'calculated by the prerun is zero.'
1253                CALL message( 'init_3d_model', 'PA0318', 1, 2, 0, 6, 0 )
1254             ENDIF
1255
1256          ENDIF
1257
1258          IF ( inflow_damping_width == 9999999.9_wp )  THEN
1259!
1260!--          Default for the transition range: one tenth of the undamped
1261!--          layer
1262             inflow_damping_width = 0.1_wp * inflow_damping_height
1263
1264          ENDIF
1265
1266          ALLOCATE( inflow_damping_factor(nzb:nzt+1) )
1267
1268          DO  k = nzb, nzt+1
1269
1270             IF ( zu(k) <= inflow_damping_height )  THEN
1271                inflow_damping_factor(k) = 1.0_wp
1272             ELSEIF ( zu(k) <= ( inflow_damping_height + inflow_damping_width ) )  THEN
1273                inflow_damping_factor(k) = 1.0_wp -                            &
1274                                           ( zu(k) - inflow_damping_height ) / &
1275                                           inflow_damping_width
1276             ELSE
1277                inflow_damping_factor(k) = 0.0_wp
1278             ENDIF
1279
1280          ENDDO
1281
1282       ENDIF
1283
1284!
1285!--    Inside buildings set velocities back to zero
1286       IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND.                &
1287            topography /= 'flat' )  THEN
1288!
1289!--       Inside buildings set velocities back to zero.
1290!--       Other scalars (pt, q, s, p, sa, ...) are ignored at present,
1291!--       maybe revise later.
1292          DO  i = nxlg, nxrg
1293             DO  j = nysg, nyng
1294                DO  k = nzb, nzt
1295                   u(k,j,i)     = MERGE( u(k,j,i), 0.0_wp,                     &
1296                                         BTEST( wall_flags_0(k,j,i), 1 ) )
1297                   v(k,j,i)     = MERGE( v(k,j,i), 0.0_wp,                     &
1298                                         BTEST( wall_flags_0(k,j,i), 2 ) )
1299                   w(k,j,i)     = MERGE( w(k,j,i), 0.0_wp,                     &
1300                                         BTEST( wall_flags_0(k,j,i), 3 ) )
1301                ENDDO
1302             ENDDO
1303          ENDDO
1304
1305       ENDIF
1306
1307!
1308!--    Calculate initial temperature field and other constants used in case
1309!--    of a sloping surface
1310       IF ( sloping_surface )  CALL init_slope
1311
1312!
1313!--    Initialize new time levels (only done in order to set boundary values
1314!--    including ghost points)
1315       pt_p = pt; u_p = u; v_p = v; w_p = w
1316       IF ( humidity )  THEN
1317          q_p = q
1318       ENDIF
1319       IF ( passive_scalar )  s_p  = s
1320!
1321!--    Allthough tendency arrays are set in prognostic_equations, they have
1322!--    have to be predefined here because they are used (but multiplied with 0)
1323!--    there before they are set.
1324       tpt_m = 0.0_wp; tu_m = 0.0_wp; tv_m = 0.0_wp; tw_m = 0.0_wp
1325       IF ( humidity )  THEN
1326          tq_m = 0.0_wp
1327       ENDIF
1328       IF ( passive_scalar )  ts_m  = 0.0_wp
1329
1330       IF ( debug_output )  CALL debug_message( 'initializing in case of restart / cyclic_fill', 'end' )
1331
1332    ELSE
1333!
1334!--    Actually this part of the programm should not be reached
1335       message_string = 'unknown initializing problem'
1336       CALL message( 'init_3d_model', 'PA0193', 1, 2, 0, 6, 0 )
1337    ENDIF
1338
1339
1340    IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1341!
1342!--    Initialize old timelevels needed for radiation boundary conditions
1343       IF ( bc_radiation_l )  THEN
1344          u_m_l(:,:,:) = u(:,:,1:2)
1345          v_m_l(:,:,:) = v(:,:,0:1)
1346          w_m_l(:,:,:) = w(:,:,0:1)
1347       ENDIF
1348       IF ( bc_radiation_r )  THEN
1349          u_m_r(:,:,:) = u(:,:,nx-1:nx)
1350          v_m_r(:,:,:) = v(:,:,nx-1:nx)
1351          w_m_r(:,:,:) = w(:,:,nx-1:nx)
1352       ENDIF
1353       IF ( bc_radiation_s )  THEN
1354          u_m_s(:,:,:) = u(:,0:1,:)
1355          v_m_s(:,:,:) = v(:,1:2,:)
1356          w_m_s(:,:,:) = w(:,0:1,:)
1357       ENDIF
1358       IF ( bc_radiation_n )  THEN
1359          u_m_n(:,:,:) = u(:,ny-1:ny,:)
1360          v_m_n(:,:,:) = v(:,ny-1:ny,:)
1361          w_m_n(:,:,:) = w(:,ny-1:ny,:)
1362       ENDIF
1363       
1364    ENDIF
1365
1366!
1367!-- Calculate the initial volume flow at the right and north boundary
1368    IF ( conserve_volume_flow )  THEN
1369
1370       IF ( use_prescribed_profile_data )  THEN
1371
1372          volume_flow_initial_l = 0.0_wp
1373          volume_flow_area_l    = 0.0_wp
1374
1375          IF ( nxr == nx )  THEN
1376             DO  j = nys, nyn
1377                DO  k = nzb+1, nzt
1378                   volume_flow_initial_l(1) = volume_flow_initial_l(1) +       &
1379                                              u_init(k) * dzw(k)               &
1380                                     * MERGE( 1.0_wp, 0.0_wp,                  &
1381                                              BTEST( wall_flags_0(k,j,nxr), 1 )&
1382                                            )
1383
1384                   volume_flow_area_l(1)    = volume_flow_area_l(1) + dzw(k)   &
1385                                     * MERGE( 1.0_wp, 0.0_wp,                  &
1386                                              BTEST( wall_flags_0(k,j,nxr), 1 )&
1387                                            )
1388                ENDDO
1389             ENDDO
1390          ENDIF
1391         
1392          IF ( nyn == ny )  THEN
1393             DO  i = nxl, nxr
1394                DO  k = nzb+1, nzt
1395                   volume_flow_initial_l(2) = volume_flow_initial_l(2) +       &
1396                                              v_init(k) * dzw(k)               &       
1397                                     * MERGE( 1.0_wp, 0.0_wp,                  &
1398                                              BTEST( wall_flags_0(k,nyn,i), 2 )&
1399                                            )
1400                   volume_flow_area_l(2)    = volume_flow_area_l(2) + dzw(k)   &       
1401                                     * MERGE( 1.0_wp, 0.0_wp,                  &
1402                                              BTEST( wall_flags_0(k,nyn,i), 2 )&
1403                                            )
1404                ENDDO
1405             ENDDO
1406          ENDIF
1407
1408#if defined( __parallel )
1409          CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&
1410                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
1411          CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1),      &
1412                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
1413
1414#else
1415          volume_flow_initial = volume_flow_initial_l
1416          volume_flow_area    = volume_flow_area_l
1417#endif 
1418
1419       ELSEIF ( TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
1420
1421          volume_flow_initial_l = 0.0_wp
1422          volume_flow_area_l    = 0.0_wp
1423
1424          IF ( nxr == nx )  THEN
1425             DO  j = nys, nyn
1426                DO  k = nzb+1, nzt
1427                   volume_flow_initial_l(1) = volume_flow_initial_l(1) +       &
1428                                              hom_sum(k,1,0) * dzw(k)          &
1429                                     * MERGE( 1.0_wp, 0.0_wp,                  &
1430                                              BTEST( wall_flags_0(k,j,nx), 1 ) &
1431                                            )
1432                   volume_flow_area_l(1)    = volume_flow_area_l(1) + dzw(k)   &
1433                                     * MERGE( 1.0_wp, 0.0_wp,                  &
1434                                              BTEST( wall_flags_0(k,j,nx), 1 ) &
1435                                            )
1436                ENDDO
1437             ENDDO
1438          ENDIF
1439         
1440          IF ( nyn == ny )  THEN
1441             DO  i = nxl, nxr
1442                DO  k = nzb+1, nzt
1443                   volume_flow_initial_l(2) = volume_flow_initial_l(2) +       &
1444                                              hom_sum(k,2,0) * dzw(k)          &       
1445                                     * MERGE( 1.0_wp, 0.0_wp,                  &
1446                                              BTEST( wall_flags_0(k,ny,i), 2 ) &
1447                                            )
1448                   volume_flow_area_l(2)    = volume_flow_area_l(2) + dzw(k)   &       
1449                                     * MERGE( 1.0_wp, 0.0_wp,                  &
1450                                              BTEST( wall_flags_0(k,ny,i), 2 ) &
1451                                            )
1452                ENDDO
1453             ENDDO
1454          ENDIF
1455
1456#if defined( __parallel )
1457          CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&
1458                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
1459          CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1),      &
1460                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
1461
1462#else
1463          volume_flow_initial = volume_flow_initial_l
1464          volume_flow_area    = volume_flow_area_l
1465#endif 
1466
1467       ELSEIF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1468
1469          volume_flow_initial_l = 0.0_wp
1470          volume_flow_area_l    = 0.0_wp
1471
1472          IF ( nxr == nx )  THEN
1473             DO  j = nys, nyn
1474                DO  k = nzb+1, nzt
1475                   volume_flow_initial_l(1) = volume_flow_initial_l(1) +       &
1476                                              u(k,j,nx) * dzw(k)               &
1477                                     * MERGE( 1.0_wp, 0.0_wp,                  &
1478                                              BTEST( wall_flags_0(k,j,nx), 1 ) &
1479                                            )
1480                   volume_flow_area_l(1)    = volume_flow_area_l(1) + dzw(k)   &
1481                                     * MERGE( 1.0_wp, 0.0_wp,                  &
1482                                              BTEST( wall_flags_0(k,j,nx), 1 ) &
1483                                            )
1484                ENDDO
1485             ENDDO
1486          ENDIF
1487         
1488          IF ( nyn == ny )  THEN
1489             DO  i = nxl, nxr
1490                DO  k = nzb+1, nzt
1491                   volume_flow_initial_l(2) = volume_flow_initial_l(2) +       &
1492                                              v(k,ny,i) * dzw(k)               &       
1493                                     * MERGE( 1.0_wp, 0.0_wp,                  &
1494                                              BTEST( wall_flags_0(k,ny,i), 2 ) &
1495                                            )
1496                   volume_flow_area_l(2)    = volume_flow_area_l(2) + dzw(k)   &       
1497                                     * MERGE( 1.0_wp, 0.0_wp,                  &
1498                                              BTEST( wall_flags_0(k,ny,i), 2 ) &
1499                                            )
1500                ENDDO
1501             ENDDO
1502          ENDIF
1503
1504#if defined( __parallel )
1505          CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&
1506                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
1507          CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1),      &
1508                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
1509
1510#else
1511          volume_flow_initial = volume_flow_initial_l
1512          volume_flow_area    = volume_flow_area_l
1513#endif 
1514
1515       ENDIF
1516
1517!
1518!--    In case of 'bulk_velocity' mode, volume_flow_initial is calculated
1519!--    from u|v_bulk instead
1520       IF ( TRIM( conserve_volume_flow_mode ) == 'bulk_velocity' )  THEN
1521          volume_flow_initial(1) = u_bulk * volume_flow_area(1)
1522          volume_flow_initial(2) = v_bulk * volume_flow_area(2)
1523       ENDIF
1524
1525    ENDIF
1526!
1527!-- In the following, surface properties can be further initialized with
1528!-- input from static driver file.
1529!-- At the moment this affects only default surfaces. For example,
1530!-- roughness length or sensible / latent heat fluxes can be initialized
1531!-- heterogeneously for default surfaces. Therefore, a generic routine
1532!-- from netcdf_data_input_mod is called to read a 2D array.
1533    IF ( input_pids_static )  THEN
1534!
1535!--    Allocate memory for possible static input
1536       ALLOCATE( tmp_2d%var(nys:nyn,nxl:nxr) )
1537       tmp_2d%var = 0.0_wp
1538!
1539!--    Open the static input file
1540#if defined( __netcdf )
1541       CALL open_read_file( TRIM( input_file_static ) //                    &
1542                            TRIM( coupling_char ),                          &
1543                            pids_static_id )
1544                           
1545       CALL inquire_num_variables( pids_static_id, num_vars_pids_static )
1546!
1547!--    Allocate memory to store variable names and read them
1548       ALLOCATE( vars_pids_static(1:num_vars_pids_static) )
1549       CALL inquire_variable_names( pids_static_id, vars_pids_static )
1550!
1551!--    Input roughness length.
1552       IF ( check_existence( vars_pids_static, 'z0' ) )  THEN
1553!
1554!--       Read _FillValue attribute
1555          CALL get_attribute( pids_static_id, char_fill, tmp_2d%fill,          &
1556                              .FALSE., 'z0' )                                 
1557!                                                                             
1558!--       Read variable                                                       
1559          CALL get_variable( pids_static_id, 'z0', tmp_2d%var,                 &
1560                             nxl, nxr, nys, nyn )                             
1561!                                                                             
1562!--       Initialize roughness length. Note, z0 will be only initialized at   
1563!--       default-type surfaces. At natural or urban z0 is implicitly         
1564!--       initialized bythe respective parameter lists.                       
1565!--       Initialize horizontal surface elements.                             
1566          CALL init_single_surface_properties( surf_def_h(0)%z0,               &
1567                                               tmp_2d%var,                     &
1568                                               surf_def_h(0)%ns,               &
1569                                               tmp_2d%fill,                    &
1570                                               surf_def_h(0)%i,                &
1571                                               surf_def_h(0)%j )               
1572!                                                                             
1573!--       Initialize roughness also at vertical surface elements.             
1574!--       Note, the actual 2D input arrays are only defined on the             
1575!--       subdomain. Therefore, pass the index arrays with their respective   
1576!--       offset values.                                                       
1577          DO  l = 0, 3                                                         
1578             CALL init_single_surface_properties(                              &
1579                                         surf_def_v(l)%z0,                     &
1580                                         tmp_2d%var,                           &
1581                                         surf_def_v(l)%ns,                     &
1582                                         tmp_2d%fill,                          &
1583                                         surf_def_v(l)%i + surf_def_v(l)%ioff, &
1584                                         surf_def_v(l)%j + surf_def_v(l)%joff )
1585          ENDDO
1586         
1587       ENDIF
1588!
1589!--    Additional variables, e.g. shf, qsws, etc, can be initialized the
1590!--    same way.
1591                           
1592!
1593!--    Finally, close the input file.
1594       CALL close_input_file( pids_static_id )
1595#endif
1596       DEALLOCATE( tmp_2d%var )
1597    ENDIF
1598!
1599!-- Finally, if random_heatflux is set, disturb shf at horizontal
1600!-- surfaces. Actually, this should be done in surface_mod, where all other
1601!-- initializations of surface quantities are done. However, this
1602!-- would create a ring dependency, hence, it is done here. Maybe delete
1603!-- disturb_heatflux and tranfer the respective code directly into the
1604!-- initialization in surface_mod.         
1605    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
1606         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
1607 
1608       IF ( use_surface_fluxes  .AND.  constant_heatflux  .AND.                &
1609            random_heatflux )  THEN
1610          IF ( surf_def_h(0)%ns >= 1 )  CALL disturb_heatflux( surf_def_h(0) )
1611          IF ( surf_lsm_h%ns    >= 1 )  CALL disturb_heatflux( surf_lsm_h    )
1612          IF ( surf_usm_h%ns    >= 1 )  CALL disturb_heatflux( surf_usm_h    )
1613       ENDIF
1614    ENDIF
1615
1616!
1617!-- Compute total sum of grid points and the mean surface level height for each
1618!-- statistic region. These are mainly used for horizontal averaging of
1619!-- turbulence statistics.
1620!-- ngp_2dh: number of grid points of a horizontal cross section through the
1621!--          respective statistic region
1622!-- ngp_3d:  number of grid points of the respective statistic region
1623    ngp_2dh_outer_l   = 0
1624    ngp_2dh_outer     = 0
1625    ngp_2dh_s_inner_l = 0
1626    ngp_2dh_s_inner   = 0
1627    ngp_2dh_l         = 0
1628    ngp_2dh           = 0
1629    ngp_3d_inner_l    = 0.0_wp
1630    ngp_3d_inner      = 0
1631    ngp_3d            = 0
1632    ngp_sums          = ( nz + 2 ) * ( pr_palm + max_pr_user )
1633
1634    mean_surface_level_height   = 0.0_wp
1635    mean_surface_level_height_l = 0.0_wp
1636!
1637!-- To do: New concept for these non-topography grid points!
1638    DO  sr = 0, statistic_regions
1639       DO  i = nxl, nxr
1640          DO  j = nys, nyn
1641             IF ( rmask(j,i,sr) == 1.0_wp )  THEN
1642!
1643!--             All xy-grid points
1644                ngp_2dh_l(sr) = ngp_2dh_l(sr) + 1
1645!
1646!--             Determine mean surface-level height. In case of downward-
1647!--             facing walls are present, more than one surface level exist.
1648!--             In this case, use the lowest surface-level height.
1649                IF ( surf_def_h(0)%start_index(j,i) <=                         &
1650                     surf_def_h(0)%end_index(j,i) )  THEN
1651                   m = surf_def_h(0)%start_index(j,i)
1652                   k = surf_def_h(0)%k(m)
1653                   mean_surface_level_height_l(sr) =                           &
1654                                       mean_surface_level_height_l(sr) + zw(k-1)
1655                ENDIF
1656                IF ( surf_lsm_h%start_index(j,i) <=                            &
1657                     surf_lsm_h%end_index(j,i) )  THEN
1658                   m = surf_lsm_h%start_index(j,i)
1659                   k = surf_lsm_h%k(m)
1660                   mean_surface_level_height_l(sr) =                           &
1661                                       mean_surface_level_height_l(sr) + zw(k-1)
1662                ENDIF
1663                IF ( surf_usm_h%start_index(j,i) <=                            &
1664                     surf_usm_h%end_index(j,i) )  THEN
1665                   m = surf_usm_h%start_index(j,i)
1666                   k = surf_usm_h%k(m)
1667                   mean_surface_level_height_l(sr) =                           &
1668                                       mean_surface_level_height_l(sr) + zw(k-1)
1669                ENDIF
1670
1671                k_surf = k - 1
1672
1673                DO  k = nzb, nzt+1
1674!
1675!--                xy-grid points above topography
1676                   ngp_2dh_outer_l(k,sr) = ngp_2dh_outer_l(k,sr)     +         &
1677                                  MERGE( 1, 0, BTEST( wall_flags_0(k,j,i), 24 ) )
1678
1679                   ngp_2dh_s_inner_l(k,sr) = ngp_2dh_s_inner_l(k,sr) +         &
1680                                  MERGE( 1, 0, BTEST( wall_flags_0(k,j,i), 22 ) )
1681
1682                ENDDO
1683!
1684!--             All grid points of the total domain above topography
1685                ngp_3d_inner_l(sr) = ngp_3d_inner_l(sr) + ( nz - k_surf + 2 )
1686
1687
1688
1689             ENDIF
1690          ENDDO
1691       ENDDO
1692    ENDDO
1693
1694    sr = statistic_regions + 1
1695#if defined( __parallel )
1696    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1697    CALL MPI_ALLREDUCE( ngp_2dh_l(0), ngp_2dh(0), sr, MPI_INTEGER, MPI_SUM,    &
1698                        comm2d, ierr )
1699    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1700    CALL MPI_ALLREDUCE( ngp_2dh_outer_l(0,0), ngp_2dh_outer(0,0), (nz+2)*sr,   &
1701                        MPI_INTEGER, MPI_SUM, comm2d, ierr )
1702    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1703    CALL MPI_ALLREDUCE( ngp_2dh_s_inner_l(0,0), ngp_2dh_s_inner(0,0),          &
1704                        (nz+2)*sr, MPI_INTEGER, MPI_SUM, comm2d, ierr )
1705    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1706    CALL MPI_ALLREDUCE( ngp_3d_inner_l(0), ngp_3d_inner_tmp(0), sr, MPI_REAL,  &
1707                        MPI_SUM, comm2d, ierr )
1708    ngp_3d_inner = INT( ngp_3d_inner_tmp, KIND = SELECTED_INT_KIND( 18 ) )
1709    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1710    CALL MPI_ALLREDUCE( mean_surface_level_height_l(0),                        &
1711                        mean_surface_level_height(0), sr, MPI_REAL,            &
1712                        MPI_SUM, comm2d, ierr )
1713    mean_surface_level_height = mean_surface_level_height / REAL( ngp_2dh )
1714#else
1715    ngp_2dh         = ngp_2dh_l
1716    ngp_2dh_outer   = ngp_2dh_outer_l
1717    ngp_2dh_s_inner = ngp_2dh_s_inner_l
1718    ngp_3d_inner    = INT( ngp_3d_inner_l, KIND = SELECTED_INT_KIND( 18 ) )
1719    mean_surface_level_height = mean_surface_level_height_l / REAL( ngp_2dh_l )
1720#endif
1721
1722    ngp_3d = INT ( ngp_2dh, KIND = SELECTED_INT_KIND( 18 ) ) * &
1723             INT ( (nz + 2 ), KIND = SELECTED_INT_KIND( 18 ) )
1724
1725!
1726!-- Set a lower limit of 1 in order to avoid zero divisions in flow_statistics,
1727!-- buoyancy, etc. A zero value will occur for cases where all grid points of
1728!-- the respective subdomain lie below the surface topography
1729    ngp_2dh_outer   = MAX( 1, ngp_2dh_outer(:,:)   ) 
1730    ngp_3d_inner    = MAX( INT(1, KIND = SELECTED_INT_KIND( 18 )),             &
1731                           ngp_3d_inner(:) )
1732    ngp_2dh_s_inner = MAX( 1, ngp_2dh_s_inner(:,:) ) 
1733
1734    DEALLOCATE( mean_surface_level_height_l, ngp_2dh_l, ngp_2dh_outer_l,       &
1735                ngp_3d_inner_l, ngp_3d_inner_tmp )
1736!
1737!-- Initialize surface forcing corresponding to large-scale forcing. Therein,
1738!-- initialize heat-fluxes, etc. via datatype. Revise it later!
1739    IF ( large_scale_forcing .AND. lsf_surf )  THEN
1740       IF ( use_surface_fluxes  .AND.  constant_heatflux )  THEN
1741          CALL ls_forcing_surf ( simulated_time )
1742       ENDIF
1743    ENDIF
1744!
1745!-- Initializae 3D offline nesting in COSMO model and read data from
1746!-- external NetCDF file.
1747    IF ( nesting_offline )  CALL nesting_offl_init
1748!
1749!-- Initialize quantities for special advections schemes
1750    CALL init_advec
1751
1752!
1753!-- Impose random perturbation on the horizontal velocity field and then
1754!-- remove the divergences from the velocity field at the initial stage
1755    IF ( create_disturbances  .AND.  disturbance_energy_limit /= 0.0_wp  .AND. &
1756         TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
1757         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
1758
1759       IF ( debug_output )  CALL debug_message( 'creating disturbances + applying pressure solver', 'start' )
1760!
1761!--    Needed for both disturb_field and pres
1762!$ACC DATA &
1763!$ACC CREATE(tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
1764!$ACC COPY(u(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
1765!$ACC COPY(v(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
1766
1767       CALL disturb_field( 'u', tend, u )
1768       CALL disturb_field( 'v', tend, v )
1769
1770!$ACC DATA &
1771!$ACC CREATE(d(nzb+1:nzt,nys:nyn,nxl:nxr)) &
1772!$ACC COPY(w(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
1773!$ACC COPY(p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
1774!$ACC COPYIN(rho_air(nzb:nzt+1), rho_air_zw(nzb:nzt+1)) &
1775!$ACC COPYIN(ddzu(1:nzt+1), ddzw(1:nzt+1)) &
1776!$ACC COPYIN(wall_flags_0(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
1777!$ACC COPYIN(bc_h(0:1)) &
1778!$ACC COPYIN(bc_h(0)%i(1:bc_h(0)%ns)) &
1779!$ACC COPYIN(bc_h(0)%j(1:bc_h(0)%ns)) &
1780!$ACC COPYIN(bc_h(0)%k(1:bc_h(0)%ns)) &
1781!$ACC COPYIN(bc_h(1)%i(1:bc_h(1)%ns)) &
1782!$ACC COPYIN(bc_h(1)%j(1:bc_h(1)%ns)) &
1783!$ACC COPYIN(bc_h(1)%k(1:bc_h(1)%ns))
1784
1785       n_sor = nsor_ini
1786       CALL pres
1787       n_sor = nsor
1788
1789!$ACC END DATA
1790!$ACC END DATA
1791
1792       IF ( debug_output )  CALL debug_message( 'creating disturbances + applying pressure solver', 'end' )
1793
1794    ENDIF
1795
1796    IF ( .NOT. ocean_mode )  THEN
1797
1798       ALLOCATE( hyp(nzb:nzt+1) )
1799       ALLOCATE( d_exner(nzb:nzt+1) )
1800       ALLOCATE( exner(nzb:nzt+1) )
1801       ALLOCATE( hyrho(nzb:nzt+1) )
1802!
1803!--    Check temperature in case of too large domain height
1804       DO  k = nzb, nzt+1
1805          IF ( ( pt_surface * exner_function(surface_pressure * 100.0_wp) - g/c_p * zu(k) ) < 0.0_wp )  THEN
1806             WRITE( message_string, * )  'absolute temperature < 0.0 at zu(', k, &
1807                                         ') = ', zu(k)
1808             CALL message( 'init_3d_model', 'PA0142', 1, 2, 0, 6, 0 )
1809          ENDIF
1810       ENDDO
1811
1812!
1813!--    Calculate vertical profile of the hydrostatic pressure (hyp)
1814       hyp    = barometric_formula(zu, pt_surface * exner_function(surface_pressure * 100.0_wp), surface_pressure * 100.0_wp)
1815       d_exner = exner_function_invers(hyp)
1816       exner = 1.0_wp / exner_function_invers(hyp)
1817       hyrho  = ideal_gas_law_rho_pt(hyp, pt_init)
1818!
1819!--    Compute reference density
1820       rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp, pt_surface * exner_function(surface_pressure * 100.0_wp))
1821
1822    ENDIF
1823
1824!
1825!-- If required, initialize particles
1826    IF ( agents_active )  CALL mas_init
1827!
1828!-- Initialization of synthetic turbulence generator
1829    IF ( use_syn_turb_gen )  CALL stg_init
1830!
1831!-- Initializing actions for all other modules
1832    CALL module_interface_init
1833!
1834!-- Initialize surface layer (done after LSM as roughness length are required
1835!-- for initialization
1836    IF ( constant_flux_layer )  CALL init_surface_layer_fluxes
1837!
1838!-- Initialize surface data output
1839    IF ( surface_output )  CALL surface_data_output_init
1840!
1841!-- Initialize the ws-scheme.   
1842    IF ( ws_scheme_sca .OR. ws_scheme_mom )  CALL ws_init
1843!
1844!-- Perform post-initializing checks for all other modules
1845    CALL module_interface_init_checks
1846
1847!
1848!-- Setting weighting factors for calculation of perturbation pressure
1849!-- and turbulent quantities from the RK substeps
1850    IF ( TRIM(timestep_scheme) == 'runge-kutta-3' )  THEN      ! for RK3-method
1851
1852       weight_substep(1) = 1._wp/6._wp
1853       weight_substep(2) = 3._wp/10._wp
1854       weight_substep(3) = 8._wp/15._wp
1855
1856       weight_pres(1)    = 1._wp/3._wp
1857       weight_pres(2)    = 5._wp/12._wp
1858       weight_pres(3)    = 1._wp/4._wp
1859
1860    ELSEIF ( TRIM(timestep_scheme) == 'runge-kutta-2' )  THEN  ! for RK2-method
1861
1862       weight_substep(1) = 1._wp/2._wp
1863       weight_substep(2) = 1._wp/2._wp
1864         
1865       weight_pres(1)    = 1._wp/2._wp
1866       weight_pres(2)    = 1._wp/2._wp       
1867
1868    ELSE                                     ! for Euler-method
1869
1870       weight_substep(1) = 1.0_wp     
1871       weight_pres(1)    = 1.0_wp                   
1872
1873    ENDIF
1874
1875!
1876!-- Initialize Rayleigh damping factors
1877    rdf    = 0.0_wp
1878    rdf_sc = 0.0_wp
1879    IF ( rayleigh_damping_factor /= 0.0_wp )  THEN
1880
1881       IF (  .NOT.  ocean_mode )  THEN
1882          DO  k = nzb+1, nzt
1883             IF ( zu(k) >= rayleigh_damping_height )  THEN
1884                rdf(k) = rayleigh_damping_factor *                             &
1885                      ( SIN( pi * 0.5_wp * ( zu(k) - rayleigh_damping_height ) &
1886                             / ( zu(nzt) - rayleigh_damping_height ) )         &
1887                      )**2
1888             ENDIF
1889          ENDDO
1890       ELSE
1891!
1892!--       In ocean mode, rayleigh damping is applied in the lower part of the
1893!--       model domain
1894          DO  k = nzt, nzb+1, -1
1895             IF ( zu(k) <= rayleigh_damping_height )  THEN
1896                rdf(k) = rayleigh_damping_factor *                             &
1897                      ( SIN( pi * 0.5_wp * ( rayleigh_damping_height - zu(k) ) &
1898                             / ( rayleigh_damping_height - zu(nzb+1) ) )       &
1899                      )**2
1900             ENDIF
1901          ENDDO
1902       ENDIF
1903
1904    ENDIF
1905    IF ( scalar_rayleigh_damping )  rdf_sc = rdf
1906
1907!
1908!-- Initialize the starting level and the vertical smoothing factor used for
1909!-- the external pressure gradient
1910    dp_smooth_factor = 1.0_wp
1911    IF ( dp_external )  THEN
1912!
1913!--    Set the starting level dp_level_ind_b only if it has not been set before
1914!--    (e.g. in init_grid).
1915       IF ( dp_level_ind_b == 0 )  THEN
1916          ind_array = MINLOC( ABS( dp_level_b - zu ) )
1917          dp_level_ind_b = ind_array(1) - 1 + nzb
1918                                        ! MINLOC uses lower array bound 1
1919       ENDIF
1920       IF ( dp_smooth )  THEN
1921          dp_smooth_factor(:dp_level_ind_b) = 0.0_wp
1922          DO  k = dp_level_ind_b+1, nzt
1923             dp_smooth_factor(k) = 0.5_wp * ( 1.0_wp + SIN( pi *               &
1924                        ( REAL( k - dp_level_ind_b, KIND=wp ) /                &
1925                          REAL( nzt - dp_level_ind_b, KIND=wp ) - 0.5_wp ) ) )
1926          ENDDO
1927       ENDIF
1928    ENDIF
1929
1930!
1931!-- Initialize damping zone for the potential temperature in case of
1932!-- non-cyclic lateral boundaries. The damping zone has the maximum value
1933!-- at the inflow boundary and decreases to zero at pt_damping_width.
1934    ptdf_x = 0.0_wp
1935    ptdf_y = 0.0_wp
1936    IF ( bc_lr_dirrad )  THEN
1937       DO  i = nxl, nxr
1938          IF ( ( i * dx ) < pt_damping_width )  THEN
1939             ptdf_x(i) = pt_damping_factor * ( SIN( pi * 0.5_wp *              &
1940                            REAL( pt_damping_width - i * dx, KIND=wp ) / (     &
1941                            REAL( pt_damping_width, KIND=wp ) ) ) )**2 
1942          ENDIF
1943       ENDDO
1944    ELSEIF ( bc_lr_raddir )  THEN
1945       DO  i = nxl, nxr
1946          IF ( ( i * dx ) > ( nx * dx - pt_damping_width ) )  THEN
1947             ptdf_x(i) = pt_damping_factor *                                   &
1948                         SIN( pi * 0.5_wp *                                    &
1949                                 ( ( i - nx ) * dx + pt_damping_width ) /      &
1950                                 REAL( pt_damping_width, KIND=wp ) )**2
1951          ENDIF
1952       ENDDO
1953    ELSEIF ( bc_ns_dirrad )  THEN
1954       DO  j = nys, nyn
1955          IF ( ( j * dy ) > ( ny * dy - pt_damping_width ) )  THEN
1956             ptdf_y(j) = pt_damping_factor *                                   &
1957                         SIN( pi * 0.5_wp *                                    &
1958                                 ( ( j - ny ) * dy + pt_damping_width ) /      &
1959                                 REAL( pt_damping_width, KIND=wp ) )**2
1960          ENDIF
1961       ENDDO
1962    ELSEIF ( bc_ns_raddir )  THEN
1963       DO  j = nys, nyn
1964          IF ( ( j * dy ) < pt_damping_width )  THEN
1965             ptdf_y(j) = pt_damping_factor *                                   &
1966                         SIN( pi * 0.5_wp *                                    &
1967                                ( pt_damping_width - j * dy ) /                &
1968                                REAL( pt_damping_width, KIND=wp ) )**2
1969          ENDIF
1970       ENDDO
1971    ENDIF
1972
1973!
1974!-- Input binary data file is not needed anymore. This line must be placed
1975!-- after call of user_init!
1976    CALL close_file( 13 )
1977!
1978!-- In case of nesting, put an barrier to assure that all parent and child
1979!-- domains finished initialization.
1980#if defined( __parallel )
1981    IF ( nested_run )  CALL MPI_BARRIER( MPI_COMM_WORLD, ierr )
1982#endif
1983
1984
1985    CALL location_message( 'model initialization', 'finished' )
1986
1987 END SUBROUTINE init_3d_model
Note: See TracBrowser for help on using the repository browser.