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

Last change on this file since 4186 was 4186, checked in by suehring, 5 years ago

Enable limitation of Obukhov length so that it does not become zero; to read input data from netcdf in init_3d_model use provided module variables instead of defining local ones; tests updated because changes in Obukhov lengths causes small differences during the initial phase of the run

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