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

Last change on this file since 4343 was 4329, checked in by motisi, 4 years ago

Renamed wall_flags_0 to wall_flags_static_0

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