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

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

Bugfix in output of time-averaged plant-canopy quanities; Output of plant-canopy data only where tall canopy is defined; land-surface model: fix wrong location strings; tests: update urban test case; all source code files: copyright update

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