source: palm/trunk/SOURCE/nesting_offl_mod.f90 @ 4270

Last change on this file since 4270 was 4270, checked in by monakurppa, 4 years ago

Updates for the aerosol module salsa

  • Reformat the emission data for the pre-processed mode (lod=2) to limit the size of the input file. Now, the emission is given as total number emission per emission category (ncat) and the number size distribution is given in variable emission_number_fracs(ncat,nbins).
  • Implement offline nesting for salsa. Additionally, offline nesting can be turn off specifically for salsa using the switch nesting_offline_salsa.
  • Remove init_aerosol_type and init_gases_type from salsa_parin and define them based on the initializing_actions
  • parameter definition removed from "season" and "season_z01" is added to parin
  • bugfix in application of index_hh after implementing the new palm_date_time_mod
  • Property svn:keywords set to Id
File size: 111.3 KB
Line 
1!> @file nesting_offl_mod.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2019 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: nesting_offl_mod.f90 4270 2019-10-23 10:46:20Z monakurppa $
27! Implement offline nesting for salsa variables.
28!
29! 4231 2019-09-12 11:22:00Z suehring
30! Bugfix in array deallocation
31!
32! 4230 2019-09-11 13:58:14Z suehring
33! Update mean chemistry profiles. These are also used for rayleigh damping.
34!
35! 4227 2019-09-10 18:04:34Z gronemeier
36! implement new palm_date_time_mod
37!
38! - Data input moved into nesting_offl_mod
39! - check rephrased
40! - time variable is now relative to time_utc_init
41! - Define module specific data type for offline nesting in nesting_offl_mod
42!
43! 4182 2019-08-22 15:20:23Z scharf
44! Corrected "Former revisions" section
45!
46! 4169 2019-08-19 13:54:35Z suehring
47! Additional check added.
48!
49! 4168 2019-08-16 13:50:17Z suehring
50! Replace function get_topography_top_index by topo_top_ind
51!
52! 4125 2019-07-29 13:31:44Z suehring
53! In order to enable netcdf parallel access, allocate dummy arrays for the
54! lateral boundary data on cores that actually do not belong to these
55! boundaries.
56!
57! 4079 2019-07-09 18:04:41Z suehring
58! - Set boundary condition for w at nzt+1 at the lateral boundaries, even
59!   though these won't enter the numerical solution. However, due to the mass
60!   conservation these values might some up to very large values which will
61!   occur in the run-control file
62! - Bugfix in offline nesting of chemical species
63! - Do not set Neumann conditions for TKE and passive scalar
64!
65! 4022 2019-06-12 11:52:39Z suehring
66! Detection of boundary-layer depth in stable boundary layer on basis of
67! boundary data improved
68! Routine for boundary-layer depth calculation renamed and made public
69!
70! 3987 2019-05-22 09:52:13Z kanani
71! Introduce alternative switch for debug output during timestepping
72!
73! 3964 2019-05-09 09:48:32Z suehring
74! Ensure that veloctiy term in calculation of bulk Richardson number does not
75! become zero
76!
77! 3937 2019-04-29 15:09:07Z suehring
78! Set boundary conditon on upper-left and upper-south grid point for the u- and
79! v-component, respectively.
80!
81! 3891 2019-04-12 17:52:01Z suehring
82! Bugfix, do not overwrite lateral and top boundary data in case of restart
83! runs.
84!
85! 3885 2019-04-11 11:29:34Z kanani
86! Changes related to global restructuring of location messages and introduction
87! of additional debug messages
88!
89!
90! Do local data exchange for chemistry variables only when boundary data is 
91! coming from dynamic file
92!
93! 3737 2019-02-12 16:57:06Z suehring
94! Introduce mesoscale nesting for chemical species
95!
96! 3705 2019-01-29 19:56:39Z suehring
97! Formatting adjustments
98!
99! 3704 2019-01-29 19:51:41Z suehring
100! Check implemented for offline nesting in child domain
101!
102! Initial Revision:
103! - separate offline nesting from large_scale_nudging_mod
104! - revise offline nesting, adjust for usage of synthetic turbulence generator
105! - adjust Rayleigh damping depending on the time-depending atmospheric
106!   conditions
107!
108!
109! Description:
110! ------------
111!> Offline nesting in larger-scale models. Boundary conditions for the simulation
112!> are read from NetCDF file and are prescribed onto the respective arrays.
113!> Further, a mass-flux correction is performed to maintain the mass balance.
114!--------------------------------------------------------------------------------!
115 MODULE nesting_offl_mod
116
117    USE arrays_3d,                                                             &
118        ONLY:  dzw,                                                            &
119               e,                                                              &
120               diss,                                                           &
121               pt,                                                             &
122               pt_init,                                                        &
123               q,                                                              &
124               q_init,                                                         &
125               rdf,                                                            &
126               rdf_sc,                                                         &
127               s,                                                              &
128               u,                                                              &
129               u_init,                                                         &
130               ug,                                                             &
131               v,                                                              &
132               v_init,                                                         &
133               vg,                                                             &
134               w,                                                              &
135               zu,                                                             &
136               zw
137
138    USE basic_constants_and_equations_mod,                                     &
139           ONLY:  g,                                                           &
140                  pi
141                 
142    USE chem_modules,                                                          &
143        ONLY:  chem_species
144
145    USE control_parameters,                                                    &
146        ONLY:  air_chemistry,                                                  & 
147               bc_dirichlet_l,                                                 & 
148               bc_dirichlet_n,                                                 &
149               bc_dirichlet_r,                                                 &
150               bc_dirichlet_s,                                                 &
151               coupling_char,                                                  &
152               dt_3d,                                                          &
153               dz,                                                             &
154               constant_diffusion,                                             &
155               child_domain,                                                   &
156               debug_output_timestep,                                          &
157               end_time,                                                       &
158               humidity,                                                       &
159               initializing_actions,                                           &
160               message_string,                                                 &
161               nesting_offline,                                                &
162               neutral,                                                        &
163               passive_scalar,                                                 &
164               rans_mode,                                                      &
165               rans_tke_e,                                                     &
166               rayleigh_damping_factor,                                        & 
167               rayleigh_damping_height,                                        &
168               salsa,                                                          &
169               spinup_time,                                                    &
170               time_since_reference_point,                                     &
171               volume_flow
172               
173    USE cpulog,                                                                &
174        ONLY:  cpu_log,                                                        &
175               log_point,                                                      &
176               log_point_s
177
178    USE grid_variables
179
180    USE indices,                                                               &
181        ONLY:  nbgp, nx, nxl, nxlg, nxlu, nxr, nxrg, ny, nys,                  &
182               nysv, nysg, nyn, nyng, nzb, nz, nzt,                            &
183               topo_top_ind,                                                   &
184               wall_flags_0
185
186    USE kinds
187
188    USE netcdf_data_input_mod,                                                 &
189        ONLY:  check_existence,                                                &
190               close_input_file,                                               &
191               get_dimension_length,                                           &
192               get_variable,                                                   &
193               get_variable_pr,                                                &
194               input_pids_dynamic,                                             &
195               inquire_num_variables,                                          &
196               inquire_variable_names,                                         &
197               input_file_dynamic,                                             &
198               num_var_pids,                                                   &
199               open_read_file,                                                 &
200               pids_id
201
202    USE palm_date_time_mod,                                                    &
203        ONLY:  get_date_time
204
205    USE pegrid
206
207    USE salsa_mod,                                                             &
208        ONLY:  salsa_nesting_offl_bc,                                          &
209               salsa_nesting_offl_init,                                        &
210               salsa_nesting_offl_input
211
212    IMPLICIT NONE
213
214!
215!-- Define data type for nesting in larger-scale models like COSMO.
216!-- Data type comprises u, v, w, pt, and q at lateral and top boundaries.
217    TYPE nest_offl_type
218
219       CHARACTER(LEN=16) ::  char_l = 'ls_forcing_left_'  !< leading substring for variables at left boundary
220       CHARACTER(LEN=17) ::  char_n = 'ls_forcing_north_' !< leading substring for variables at north boundary 
221       CHARACTER(LEN=17) ::  char_r = 'ls_forcing_right_' !< leading substring for variables at right boundary 
222       CHARACTER(LEN=17) ::  char_s = 'ls_forcing_south_' !< leading substring for variables at south boundary
223       CHARACTER(LEN=15) ::  char_t = 'ls_forcing_top_'   !< leading substring for variables at top boundary
224
225       CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE ::  var_names         !< list of variable in dynamic input file
226       CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE ::  var_names_chem_l  !< names of mesoscale nested chemistry variables at left boundary
227       CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE ::  var_names_chem_n  !< names of mesoscale nested chemistry variables at north boundary
228       CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE ::  var_names_chem_r  !< names of mesoscale nested chemistry variables at right boundary
229       CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE ::  var_names_chem_s  !< names of mesoscale nested chemistry variables at south boundary
230       CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE ::  var_names_chem_t  !< names of mesoscale nested chemistry variables at top boundary
231
232       INTEGER(iwp) ::  nt     !< number of time levels in dynamic input file
233       INTEGER(iwp) ::  nzu    !< number of vertical levels on scalar grid in dynamic input file
234       INTEGER(iwp) ::  nzw    !< number of vertical levels on w grid in dynamic input file
235       INTEGER(iwp) ::  tind   !< time index for reference time in mesoscale-offline nesting
236       INTEGER(iwp) ::  tind_p !< time index for following time in mesoscale-offline nesting
237
238       LOGICAL      ::  init         = .FALSE. !< flag indicating that offline nesting is already initialized
239
240       LOGICAL, DIMENSION(:), ALLOCATABLE ::  chem_from_file_l !< flags inidicating whether left boundary data for chemistry is in dynamic input file 
241       LOGICAL, DIMENSION(:), ALLOCATABLE ::  chem_from_file_n !< flags inidicating whether north boundary data for chemistry is in dynamic input file
242       LOGICAL, DIMENSION(:), ALLOCATABLE ::  chem_from_file_r !< flags inidicating whether right boundary data for chemistry is in dynamic input file
243       LOGICAL, DIMENSION(:), ALLOCATABLE ::  chem_from_file_s !< flags inidicating whether south boundary data for chemistry is in dynamic input file
244       LOGICAL, DIMENSION(:), ALLOCATABLE ::  chem_from_file_t !< flags inidicating whether top boundary data for chemistry is in dynamic input file
245
246       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surface_pressure !< time dependent surface pressure
247       REAL(wp), DIMENSION(:), ALLOCATABLE ::  time             !< time levels in dynamic input file
248       REAL(wp), DIMENSION(:), ALLOCATABLE ::  zu_atmos         !< vertical levels at scalar grid in dynamic input file
249       REAL(wp), DIMENSION(:), ALLOCATABLE ::  zw_atmos         !< vertical levels at w grid in dynamic input file
250
251       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ug         !< domain-averaged geostrophic component
252       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  vg         !< domain-averaged geostrophic component
253
254       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  u_left   !< u-component at left boundary
255       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  v_left   !< v-component at left boundary
256       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  w_left   !< w-component at left boundary
257       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  q_left   !< mixing ratio at left boundary
258       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pt_left  !< potentital temperautre at left boundary
259
260       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  u_north  !< u-component at north boundary
261       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  v_north  !< v-component at north boundary
262       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  w_north  !< w-component at north boundary
263       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  q_north  !< mixing ratio at north boundary
264       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pt_north !< potentital temperautre at north boundary
265
266       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  u_right  !< u-component at right boundary
267       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  v_right  !< v-component at right boundary
268       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  w_right  !< w-component at right boundary
269       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  q_right  !< mixing ratio at right boundary
270       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pt_right !< potentital temperautre at right boundary
271
272       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  u_south  !< u-component at south boundary
273       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  v_south  !< v-component at south boundary
274       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  w_south  !< w-component at south boundary
275       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  q_south  !< mixing ratio at south boundary
276       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pt_south !< potentital temperautre at south boundary
277
278       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  u_top    !< u-component at top boundary
279       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  v_top    !< v-component at top boundary
280       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  w_top    !< w-component at top boundary
281       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  q_top    !< mixing ratio at top boundary
282       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pt_top   !< potentital temperautre at top boundary
283       
284       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  chem_left   !< chemical species at left boundary
285       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  chem_north  !< chemical species at left boundary
286       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  chem_right  !< chemical species at left boundary
287       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  chem_south  !< chemical species at left boundary
288       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  chem_top    !< chemical species at left boundary
289
290    END TYPE nest_offl_type
291
292    REAL(wp) ::  fac_dt              !< interpolation factor
293    REAL(wp) ::  time_utc_init       !< time in seconds-of-day of origin_date_time
294    REAL(wp) ::  zi_ribulk = 0.0_wp  !< boundary-layer depth according to bulk Richardson criterion, i.e. the height where Ri_bulk exceeds the critical
295                                     !< bulk Richardson number of 0.2
296
297    TYPE(nest_offl_type) ::  nest_offl  !< data structure for data input at lateral and top boundaries (provided by Inifor)
298   
299    SAVE
300    PRIVATE
301!
302!-- Public subroutines
303    PUBLIC nesting_offl_bc,                                                    &
304           nesting_offl_calc_zi,                                               &
305           nesting_offl_check_parameters,                                      &
306           nesting_offl_geostrophic_wind,                                      &
307           nesting_offl_header,                                                &
308           nesting_offl_init,                                                  &
309           nesting_offl_input,                                                 &
310           nesting_offl_interpolation_factor,                                  &
311           nesting_offl_mass_conservation,                                     &
312           nesting_offl_parin 
313!
314!-- Public variables
315    PUBLIC zi_ribulk   
316
317    INTERFACE nesting_offl_bc
318       MODULE PROCEDURE nesting_offl_bc
319    END INTERFACE nesting_offl_bc
320   
321    INTERFACE nesting_offl_calc_zi
322       MODULE PROCEDURE nesting_offl_calc_zi
323    END INTERFACE nesting_offl_calc_zi
324   
325    INTERFACE nesting_offl_check_parameters
326       MODULE PROCEDURE nesting_offl_check_parameters
327    END INTERFACE nesting_offl_check_parameters
328
329    INTERFACE nesting_offl_geostrophic_wind
330       MODULE PROCEDURE nesting_offl_geostrophic_wind
331    END INTERFACE nesting_offl_geostrophic_wind
332   
333    INTERFACE nesting_offl_header
334       MODULE PROCEDURE nesting_offl_header
335    END INTERFACE nesting_offl_header
336   
337    INTERFACE nesting_offl_init
338       MODULE PROCEDURE nesting_offl_init
339    END INTERFACE nesting_offl_init
340
341    INTERFACE nesting_offl_input
342       MODULE PROCEDURE nesting_offl_input
343    END INTERFACE nesting_offl_input
344
345    INTERFACE nesting_offl_interpolation_factor
346       MODULE PROCEDURE nesting_offl_interpolation_factor
347    END INTERFACE nesting_offl_interpolation_factor
348           
349    INTERFACE nesting_offl_mass_conservation
350       MODULE PROCEDURE nesting_offl_mass_conservation
351    END INTERFACE nesting_offl_mass_conservation
352   
353    INTERFACE nesting_offl_parin
354       MODULE PROCEDURE nesting_offl_parin
355    END INTERFACE nesting_offl_parin
356
357 CONTAINS
358
359!------------------------------------------------------------------------------!
360! Description:
361! ------------
362!> Reads data at lateral and top boundaries derived from larger-scale model.
363!------------------------------------------------------------------------------!
364    SUBROUTINE nesting_offl_input
365
366       INTEGER(iwp) ::  n   !< running index for chemistry variables
367       INTEGER(iwp) ::  t   !< running index time dimension
368
369!
370!--    Initialize INIFOR forcing in first call.
371       IF ( .NOT. nest_offl%init )  THEN
372#if defined ( __netcdf )
373!
374!--       Open file in read-only mode
375          CALL open_read_file( TRIM( input_file_dynamic ) //                   &
376                               TRIM( coupling_char ), pids_id )
377!
378!--       At first, inquire all variable names.
379          CALL inquire_num_variables( pids_id, num_var_pids )
380!
381!--       Allocate memory to store variable names.
382          ALLOCATE( nest_offl%var_names(1:num_var_pids) )
383          CALL inquire_variable_names( pids_id, nest_offl%var_names )
384!
385!--       Read time dimension, allocate memory and finally read time array
386          CALL get_dimension_length( pids_id, nest_offl%nt, 'time' )
387
388          IF ( check_existence( nest_offl%var_names, 'time' ) )  THEN
389             ALLOCATE( nest_offl%time(0:nest_offl%nt-1) )
390             CALL get_variable( pids_id, 'time', nest_offl%time )
391          ENDIF
392!
393!--       Read vertical dimension of scalar und w grid
394          CALL get_dimension_length( pids_id, nest_offl%nzu, 'z' )
395          CALL get_dimension_length( pids_id, nest_offl%nzw, 'zw' )
396
397          IF ( check_existence( nest_offl%var_names, 'z' ) )  THEN
398             ALLOCATE( nest_offl%zu_atmos(1:nest_offl%nzu) )
399             CALL get_variable( pids_id, 'z', nest_offl%zu_atmos )
400          ENDIF
401          IF ( check_existence( nest_offl%var_names, 'zw' ) )  THEN
402             ALLOCATE( nest_offl%zw_atmos(1:nest_offl%nzw) )
403             CALL get_variable( pids_id, 'zw', nest_offl%zw_atmos )
404          ENDIF
405!
406!--       Read surface pressure
407          IF ( check_existence( nest_offl%var_names,                           &
408                                'surface_forcing_surface_pressure' ) )  THEN
409             ALLOCATE( nest_offl%surface_pressure(0:nest_offl%nt-1) )
410             CALL get_variable( pids_id,                                       &
411                                'surface_forcing_surface_pressure',            &
412                                nest_offl%surface_pressure )
413          ENDIF
414!
415!--       Close input file
416          CALL close_input_file( pids_id )
417#endif
418       ENDIF
419!
420!--    Check if dynamic driver data input is required.
421       IF ( nest_offl%time(nest_offl%tind_p) <=                                &
422            MAX( time_since_reference_point, 0.0_wp) + time_utc_init  .OR.     &
423            .NOT.  nest_offl%init )  THEN
424          CONTINUE
425!
426!--    Return otherwise
427       ELSE
428          RETURN
429       ENDIF
430!
431!--    CPU measurement
432       CALL cpu_log( log_point_s(86), 'NetCDF input forcing', 'start' )
433
434!
435!--    Obtain time index for current point in time. Note, the time coordinate
436!--    in the input file is relative to time_utc_init. Since time_since_...
437!--    is negativ when spinup is used, use MAX function to obtain correct
438!--    time at the beginning.
439       nest_offl%tind = MINLOC( ABS( nest_offl%time - (                        &
440                                     time_utc_init +                           &
441                                     MAX( time_since_reference_point, 0.0_wp) )&
442                                   ), DIM = 1 ) - 1
443       nest_offl%tind_p = nest_offl%tind + 1
444!
445!--    Open file in read-only mode
446#if defined ( __netcdf )
447       CALL open_read_file( TRIM( input_file_dynamic ) //                      &
448                            TRIM( coupling_char ), pids_id )
449!
450!--    Read geostrophic wind components
451       DO  t = nest_offl%tind, nest_offl%tind_p
452          CALL get_variable_pr( pids_id, 'ls_forcing_ug', t+1,                 &
453                                nest_offl%ug(t-nest_offl%tind,nzb+1:nzt) )
454          CALL get_variable_pr( pids_id, 'ls_forcing_vg', t+1,                 &
455                                nest_offl%vg(t-nest_offl%tind,nzb+1:nzt) )
456       ENDDO
457!
458!--    Read data at lateral and top boundaries. Please note, at left and
459!--    right domain boundary, yz-layers are read for u, v, w, pt and q.
460!--    For the v-component, the data starts at nysv, while for the other
461!--    quantities the data starts at nys. This is equivalent at the north
462!--    and south domain boundary for the u-component.
463!--    Note, lateral data is also accessed by parallel IO, which is the reason
464!--    why different arguments are passed depending on the boundary control
465!--    flags. Cores that do not belong to the respective boundary just make
466!--    a dummy read with count = 0, just in order to participate the collective
467!--    operation.
468!--    Read data for western boundary   
469       CALL get_variable( pids_id, 'ls_forcing_left_u',                        &
470                          nest_offl%u_left,                                    & ! array to be read
471                          MERGE( nys+1, 1, bc_dirichlet_l),                    & ! start index y direction
472                          MERGE( nzb+1, 1, bc_dirichlet_l),                    & ! start index z direction
473                          MERGE( nest_offl%tind+1, 1, bc_dirichlet_l),         & ! start index time dimension
474                          MERGE( nyn-nys+1, 0, bc_dirichlet_l),                & ! number of elements along y
475                          MERGE( nest_offl%nzu, 0, bc_dirichlet_l),            & ! number of elements alogn z
476                          MERGE( 2, 0, bc_dirichlet_l),                        & ! number of time steps (2 or 0)
477                          .TRUE. )                                               ! parallel IO when compiled accordingly
478     
479       CALL get_variable( pids_id, 'ls_forcing_left_v',                        &
480                          nest_offl%v_left,                                    &
481                          MERGE( nysv, 1, bc_dirichlet_l),                     &
482                          MERGE( nzb+1, 1, bc_dirichlet_l),                    &
483                          MERGE( nest_offl%tind+1, 1, bc_dirichlet_l),         &
484                          MERGE( nyn-nysv+1, 0, bc_dirichlet_l),               &
485                          MERGE( nest_offl%nzu, 0, bc_dirichlet_l),            &
486                          MERGE( 2, 0, bc_dirichlet_l),                        &
487                          .TRUE. )                                       
488
489       CALL get_variable( pids_id, 'ls_forcing_left_w',                        &
490                          nest_offl%w_left,                                    &
491                          MERGE( nys+1, 1, bc_dirichlet_l),                    &
492                          MERGE( nzb+1, 1, bc_dirichlet_l),                    &
493                          MERGE( nest_offl%tind+1, 1, bc_dirichlet_l),         &
494                          MERGE( nyn-nys+1, 0, bc_dirichlet_l),                &
495                          MERGE( nest_offl%nzw, 0, bc_dirichlet_l),            &
496                          MERGE( 2, 0, bc_dirichlet_l),                        &
497                          .TRUE. )   
498
499       IF ( .NOT. neutral )  THEN
500          CALL get_variable( pids_id, 'ls_forcing_left_pt',                    &
501                             nest_offl%pt_left,                                &
502                             MERGE( nys+1, 1, bc_dirichlet_l),                 &
503                             MERGE( nzb+1, 1, bc_dirichlet_l),                 &
504                             MERGE( nest_offl%tind+1, 1, bc_dirichlet_l),      &
505                             MERGE( nyn-nys+1, 0, bc_dirichlet_l),             &
506                             MERGE( nest_offl%nzu, 0, bc_dirichlet_l),         &
507                             MERGE( 2, 0, bc_dirichlet_l),                     &
508                             .TRUE. )
509       ENDIF
510
511       IF ( humidity )  THEN
512          CALL get_variable( pids_id, 'ls_forcing_left_qv',                    &
513                             nest_offl%q_left,                                 &
514                             MERGE( nys+1, 1, bc_dirichlet_l),                 &
515                             MERGE( nzb+1, 1, bc_dirichlet_l),                 &
516                             MERGE( nest_offl%tind+1, 1, bc_dirichlet_l),      &
517                             MERGE( nyn-nys+1, 0, bc_dirichlet_l),             &
518                             MERGE( nest_offl%nzu, 0, bc_dirichlet_l),         &
519                             MERGE( 2, 0, bc_dirichlet_l),                     &
520                             .TRUE. )
521       ENDIF
522       
523       IF ( air_chemistry )  THEN
524          DO  n = 1, UBOUND(nest_offl%var_names_chem_l, 1)
525             IF ( check_existence( nest_offl%var_names,                        &
526                                   nest_offl%var_names_chem_l(n) ) )  THEN 
527                CALL get_variable( pids_id,                                    &
528                           TRIM( nest_offl%var_names_chem_l(n) ),              &
529                           nest_offl%chem_left(:,:,:,n),                       &
530                           MERGE( nys+1, 1, bc_dirichlet_l),                   &
531                           MERGE( nzb+1, 1, bc_dirichlet_l),                   &
532                           MERGE( nest_offl%tind+1, 1, bc_dirichlet_l),        &
533                           MERGE( nyn-nys+1, 0, bc_dirichlet_l),               &
534                           MERGE( nest_offl%nzu, 0, bc_dirichlet_l),           &
535                           MERGE( 2, 0, bc_dirichlet_l),                       &
536                           .TRUE. )
537                nest_offl%chem_from_file_l(n) = .TRUE.
538             ENDIF
539          ENDDO
540       ENDIF
541!
542!--    Read data for eastern boundary   
543       CALL get_variable( pids_id, 'ls_forcing_right_u',                       &
544                          nest_offl%u_right,                                   &
545                          MERGE( nys+1, 1, bc_dirichlet_r),                    &
546                          MERGE( nzb+1, 1, bc_dirichlet_r),                    &
547                          MERGE( nest_offl%tind+1, 1, bc_dirichlet_r),         &
548                          MERGE( nyn-nys+1, 0, bc_dirichlet_r),                &
549                          MERGE( nest_offl%nzu, 0, bc_dirichlet_r),            &
550                          MERGE( 2, 0, bc_dirichlet_r),                        &
551                          .TRUE. )                                             
552     
553       CALL get_variable( pids_id, 'ls_forcing_right_v',                       &
554                          nest_offl%v_right,                                   &
555                          MERGE( nysv, 1, bc_dirichlet_r),                     &
556                          MERGE( nzb+1, 1, bc_dirichlet_r),                    &
557                          MERGE( nest_offl%tind+1, 1, bc_dirichlet_r),         &
558                          MERGE( nyn-nysv+1, 0, bc_dirichlet_r),               &
559                          MERGE( nest_offl%nzu, 0, bc_dirichlet_r),            &
560                          MERGE( 2, 0, bc_dirichlet_r),                        &
561                          .TRUE. )                                             
562
563       CALL get_variable( pids_id, 'ls_forcing_right_w',                       &
564                          nest_offl%w_right,                                   &
565                          MERGE( nys+1, 1, bc_dirichlet_r),                    &
566                          MERGE( nzb+1, 1, bc_dirichlet_r),                    &
567                          MERGE( nest_offl%tind+1, 1, bc_dirichlet_r),         &
568                          MERGE( nyn-nys+1, 0, bc_dirichlet_r),                &
569                          MERGE( nest_offl%nzw, 0, bc_dirichlet_r),            &
570                          MERGE( 2, 0, bc_dirichlet_r),                        &
571                          .TRUE. )   
572
573       IF ( .NOT. neutral )  THEN
574          CALL get_variable( pids_id, 'ls_forcing_right_pt',                   &
575                             nest_offl%pt_right,                               &
576                             MERGE( nys+1, 1, bc_dirichlet_r),                 &
577                             MERGE( nzb+1, 1, bc_dirichlet_r),                 &
578                             MERGE( nest_offl%tind+1, 1, bc_dirichlet_r),      &
579                             MERGE( nyn-nys+1, 0, bc_dirichlet_r),             &
580                             MERGE( nest_offl%nzu, 0, bc_dirichlet_r),         &
581                             MERGE( 2, 0, bc_dirichlet_r),                     &
582                             .TRUE. )
583       ENDIF
584
585       IF ( humidity )  THEN
586          CALL get_variable( pids_id, 'ls_forcing_right_qv',                   &
587                             nest_offl%q_right,                                &
588                             MERGE( nys+1, 1, bc_dirichlet_r),                 &
589                             MERGE( nzb+1, 1, bc_dirichlet_r),                 &
590                             MERGE( nest_offl%tind+1, 1, bc_dirichlet_r),      &
591                             MERGE( nyn-nys+1, 0, bc_dirichlet_r),             &
592                             MERGE( nest_offl%nzu, 0, bc_dirichlet_r),         &
593                             MERGE( 2, 0, bc_dirichlet_r),                     &
594                             .TRUE. )
595       ENDIF
596       
597       IF ( air_chemistry )  THEN
598          DO  n = 1, UBOUND(nest_offl%var_names_chem_r, 1)
599             IF ( check_existence( nest_offl%var_names,                        &
600                                   nest_offl%var_names_chem_r(n) ) )  THEN     
601                CALL get_variable( pids_id,                                    &
602                           TRIM( nest_offl%var_names_chem_r(n) ),              &
603                           nest_offl%chem_right(:,:,:,n),                      &
604                           MERGE( nys+1, 1, bc_dirichlet_r),                   &
605                           MERGE( nzb+1, 1, bc_dirichlet_r),                   &
606                           MERGE( nest_offl%tind+1, 1, bc_dirichlet_r),        &
607                           MERGE( nyn-nys+1, 0, bc_dirichlet_r),               &
608                           MERGE( nest_offl%nzu, 0, bc_dirichlet_r),           &
609                           MERGE( 2, 0, bc_dirichlet_r),                       &
610                           .TRUE. )
611                nest_offl%chem_from_file_r(n) = .TRUE.
612             ENDIF
613          ENDDO
614       ENDIF
615!
616!--    Read data for northern boundary
617       CALL get_variable( pids_id, 'ls_forcing_north_u',                       & ! array to be read
618                          nest_offl%u_north,                                   & ! start index x direction
619                          MERGE( nxlu, 1, bc_dirichlet_n ),                    & ! start index z direction
620                          MERGE( nzb+1, 1, bc_dirichlet_n ),                   & ! start index time dimension
621                          MERGE( nest_offl%tind+1, 1, bc_dirichlet_n ),        & ! number of elements along x
622                          MERGE( nxr-nxlu+1, 0, bc_dirichlet_n ),              & ! number of elements alogn z
623                          MERGE( nest_offl%nzu, 0, bc_dirichlet_n ),           & ! number of time steps (2 or 0)
624                          MERGE( 2, 0, bc_dirichlet_n ),                       & ! parallel IO when compiled accordingly
625                          .TRUE. )                                             
626                                                                               
627       CALL get_variable( pids_id, 'ls_forcing_north_v',                       & ! array to be read
628                          nest_offl%v_north,                                   & ! start index x direction
629                          MERGE( nxl+1, 1, bc_dirichlet_n ),                   & ! start index z direction
630                          MERGE( nzb+1, 1, bc_dirichlet_n ),                   & ! start index time dimension
631                          MERGE( nest_offl%tind+1, 1, bc_dirichlet_n ),        & ! number of elements along x
632                          MERGE( nxr-nxl+1, 0, bc_dirichlet_n ),               & ! number of elements alogn z
633                          MERGE( nest_offl%nzu, 0, bc_dirichlet_n ),           & ! number of time steps (2 or 0)
634                          MERGE( 2, 0, bc_dirichlet_n ),                       & ! parallel IO when compiled accordingly
635                          .TRUE. )                                             
636                                                                               
637       CALL get_variable( pids_id, 'ls_forcing_north_w',                       & ! array to be read
638                          nest_offl%w_north,                                   & ! start index x direction
639                          MERGE( nxl+1, 1, bc_dirichlet_n ),                   & ! start index z direction
640                          MERGE( nzb+1, 1, bc_dirichlet_n ),                   & ! start index time dimension
641                          MERGE( nest_offl%tind+1, 1, bc_dirichlet_n ),        & ! number of elements along x
642                          MERGE( nxr-nxl+1, 0, bc_dirichlet_n ),               & ! number of elements alogn z
643                          MERGE( nest_offl%nzw, 0, bc_dirichlet_n ),           & ! number of time steps (2 or 0)
644                          MERGE( 2, 0, bc_dirichlet_n ),                       & ! parallel IO when compiled accordingly
645                          .TRUE. )                                             
646                                                                               
647       IF ( .NOT. neutral )  THEN                                             
648          CALL get_variable( pids_id, 'ls_forcing_north_pt',                   & ! array to be read
649                             nest_offl%pt_north,                               & ! start index x direction
650                             MERGE( nxl+1, 1, bc_dirichlet_n ),                & ! start index z direction
651                             MERGE( nzb+1, 1, bc_dirichlet_n ),                & ! start index time dimension
652                             MERGE( nest_offl%tind+1, 1, bc_dirichlet_n ),     & ! number of elements along x
653                             MERGE( nxr-nxl+1, 0, bc_dirichlet_n ),            & ! number of elements alogn z
654                             MERGE( nest_offl%nzu, 0, bc_dirichlet_n ),        & ! number of time steps (2 or 0)
655                             MERGE( 2, 0, bc_dirichlet_n ),                    & ! parallel IO when compiled accordingly
656                             .TRUE. )                                             
657       ENDIF                                                                   
658       IF ( humidity )  THEN                                                   
659          CALL get_variable( pids_id, 'ls_forcing_north_qv',                   & ! array to be read
660                             nest_offl%q_north,                                & ! start index x direction
661                             MERGE( nxl+1, 1, bc_dirichlet_n ),                & ! start index z direction
662                             MERGE( nzb+1, 1, bc_dirichlet_n ),                & ! start index time dimension
663                             MERGE( nest_offl%tind+1, 1, bc_dirichlet_n ),     & ! number of elements along x
664                             MERGE( nxr-nxl+1, 0, bc_dirichlet_n ),            & ! number of elements alogn z
665                             MERGE( nest_offl%nzu, 0, bc_dirichlet_n ),        & ! number of time steps (2 or 0)
666                             MERGE( 2, 0, bc_dirichlet_n ),                    & ! parallel IO when compiled accordingly
667                             .TRUE. )                                             
668       ENDIF                                                                   
669                                                                               
670       IF ( air_chemistry )  THEN                                             
671          DO  n = 1, UBOUND(nest_offl%var_names_chem_n, 1)                     
672             IF ( check_existence( nest_offl%var_names,                        &
673                                   nest_offl%var_names_chem_n(n) ) )  THEN     
674                CALL get_variable( pids_id,                                    &
675                           TRIM( nest_offl%var_names_chem_n(n) ),              &
676                           nest_offl%chem_north(:,:,:,n),                      &
677                           MERGE( nxl+1, 1, bc_dirichlet_n ),                  &
678                           MERGE( nzb+1, 1, bc_dirichlet_n ),                  &
679                           MERGE( nest_offl%tind+1, 1, bc_dirichlet_n ),       &
680                           MERGE( nxr-nxl+1, 0, bc_dirichlet_n ),              &
681                           MERGE( nest_offl%nzu, 0, bc_dirichlet_n ),          &
682                           MERGE( 2, 0, bc_dirichlet_n ),                      &
683                           .TRUE. )
684                nest_offl%chem_from_file_n(n) = .TRUE.
685             ENDIF
686          ENDDO
687       ENDIF
688!
689!--    Read data for southern boundary
690       CALL get_variable( pids_id, 'ls_forcing_south_u',                       & ! array to be read
691                          nest_offl%u_south,                                   & ! start index x direction
692                          MERGE( nxlu, 1, bc_dirichlet_s ),                    & ! start index z direction
693                          MERGE( nzb+1, 1, bc_dirichlet_s ),                   & ! start index time dimension
694                          MERGE( nest_offl%tind+1, 1, bc_dirichlet_s ),        & ! number of elements along x
695                          MERGE( nxr-nxlu+1, 0, bc_dirichlet_s ),              & ! number of elements alogn z
696                          MERGE( nest_offl%nzu, 0, bc_dirichlet_s ),           & ! number of time steps (2 or 0)
697                          MERGE( 2, 0, bc_dirichlet_s ),                       & ! parallel IO when compiled accordingly
698                          .TRUE. )                                             
699                                                                               
700       CALL get_variable( pids_id, 'ls_forcing_south_v',                       & ! array to be read
701                          nest_offl%v_south,                                   & ! start index x direction
702                          MERGE( nxl+1, 1, bc_dirichlet_s ),                   & ! start index z direction
703                          MERGE( nzb+1, 1, bc_dirichlet_s ),                   & ! start index time dimension
704                          MERGE( nest_offl%tind+1, 1, bc_dirichlet_s ),        & ! number of elements along x
705                          MERGE( nxr-nxl+1, 0, bc_dirichlet_s ),               & ! number of elements alogn z
706                          MERGE( nest_offl%nzu, 0, bc_dirichlet_s ),           & ! number of time steps (2 or 0)
707                          MERGE( 2, 0, bc_dirichlet_s ),                       & ! parallel IO when compiled accordingly
708                          .TRUE. )                                             
709                                                                               
710       CALL get_variable( pids_id, 'ls_forcing_south_w',                       & ! array to be read
711                          nest_offl%w_south,                                   & ! start index x direction
712                          MERGE( nxl+1, 1, bc_dirichlet_s ),                   & ! start index z direction
713                          MERGE( nzb+1, 1, bc_dirichlet_s ),                   & ! start index time dimension
714                          MERGE( nest_offl%tind+1, 1, bc_dirichlet_s ),        & ! number of elements along x
715                          MERGE( nxr-nxl+1, 0, bc_dirichlet_s ),               & ! number of elements alogn z
716                          MERGE( nest_offl%nzw, 0, bc_dirichlet_s ),           & ! number of time steps (2 or 0)
717                          MERGE( 2, 0, bc_dirichlet_s ),                       & ! parallel IO when compiled accordingly
718                          .TRUE. )                                             
719                                                                               
720       IF ( .NOT. neutral )  THEN                                             
721          CALL get_variable( pids_id, 'ls_forcing_south_pt',                   & ! array to be read
722                             nest_offl%pt_south,                               & ! start index x direction
723                             MERGE( nxl+1, 1, bc_dirichlet_s ),                & ! start index z direction
724                             MERGE( nzb+1, 1, bc_dirichlet_s ),                & ! start index time dimension
725                             MERGE( nest_offl%tind+1, 1, bc_dirichlet_s ),     & ! number of elements along x
726                             MERGE( nxr-nxl+1, 0, bc_dirichlet_s ),            & ! number of elements alogn z
727                             MERGE( nest_offl%nzu, 0, bc_dirichlet_s ),        & ! number of time steps (2 or 0)
728                             MERGE( 2, 0, bc_dirichlet_s ),                    & ! parallel IO when compiled accordingly
729                             .TRUE. )                                             
730       ENDIF                                                                   
731       IF ( humidity )  THEN                                                   
732          CALL get_variable( pids_id, 'ls_forcing_south_qv',                   & ! array to be read
733                             nest_offl%q_south,                                & ! start index x direction
734                             MERGE( nxl+1, 1, bc_dirichlet_s ),                & ! start index z direction
735                             MERGE( nzb+1, 1, bc_dirichlet_s ),                & ! start index time dimension
736                             MERGE( nest_offl%tind+1, 1, bc_dirichlet_s ),     & ! number of elements along x
737                             MERGE( nxr-nxl+1, 0, bc_dirichlet_s ),            & ! number of elements alogn z
738                             MERGE( nest_offl%nzu, 0, bc_dirichlet_s ),        & ! number of time steps (2 or 0)
739                             MERGE( 2, 0, bc_dirichlet_s ),                    & ! parallel IO when compiled accordingly
740                             .TRUE. )                                             
741       ENDIF                                                                   
742                                                                               
743       IF ( air_chemistry )  THEN                                             
744          DO  n = 1, UBOUND(nest_offl%var_names_chem_s, 1)                     
745             IF ( check_existence( nest_offl%var_names,                        &
746                                   nest_offl%var_names_chem_s(n) ) )  THEN     
747                CALL get_variable( pids_id,                                    &
748                           TRIM( nest_offl%var_names_chem_s(n) ),              &
749                           nest_offl%chem_south(:,:,:,n),                      &
750                           MERGE( nxl+1, 1, bc_dirichlet_s ),                  &
751                           MERGE( nzb+1, 1, bc_dirichlet_s ),                  &
752                           MERGE( nest_offl%tind+1, 1, bc_dirichlet_s ),       &
753                           MERGE( nxr-nxl+1, 0, bc_dirichlet_s ),              &
754                           MERGE( nest_offl%nzu, 0, bc_dirichlet_s ),          &
755                           MERGE( 2, 0, bc_dirichlet_s ),                      &
756                           .TRUE. )
757                nest_offl%chem_from_file_s(n) = .TRUE.
758             ENDIF
759          ENDDO
760       ENDIF
761!
762!--    Top boundary
763       CALL get_variable( pids_id, 'ls_forcing_top_u',                         &
764                             nest_offl%u_top(0:1,nys:nyn,nxlu:nxr),            &
765                             nxlu, nys+1, nest_offl%tind+1,                    &
766                             nxr-nxlu+1, nyn-nys+1, 2, .TRUE. )
767
768       CALL get_variable( pids_id, 'ls_forcing_top_v',                         &
769                             nest_offl%v_top(0:1,nysv:nyn,nxl:nxr),            &
770                             nxl+1, nysv, nest_offl%tind+1,                    &
771                             nxr-nxl+1, nyn-nysv+1, 2, .TRUE. )
772                             
773       CALL get_variable( pids_id, 'ls_forcing_top_w',                         &
774                             nest_offl%w_top(0:1,nys:nyn,nxl:nxr),             &
775                             nxl+1, nys+1, nest_offl%tind+1,                   &
776                             nxr-nxl+1, nyn-nys+1, 2, .TRUE. )
777                             
778       IF ( .NOT. neutral )  THEN
779          CALL get_variable( pids_id, 'ls_forcing_top_pt',                     &
780                                nest_offl%pt_top(0:1,nys:nyn,nxl:nxr),         &
781                                nxl+1, nys+1, nest_offl%tind+1,                &
782                                nxr-nxl+1, nyn-nys+1, 2, .TRUE. )
783       ENDIF
784       IF ( humidity )  THEN
785          CALL get_variable( pids_id, 'ls_forcing_top_qv',                     &
786                                nest_offl%q_top(0:1,nys:nyn,nxl:nxr),          &
787                                nxl+1, nys+1, nest_offl%tind+1,                &
788                                nxr-nxl+1, nyn-nys+1, 2, .TRUE. )
789       ENDIF
790       
791       IF ( air_chemistry )  THEN
792          DO  n = 1, UBOUND(nest_offl%var_names_chem_t, 1)
793             IF ( check_existence( nest_offl%var_names,                        &
794                                   nest_offl%var_names_chem_t(n) ) )  THEN     
795                CALL get_variable( pids_id,                                    &
796                              TRIM( nest_offl%var_names_chem_t(n) ),           &
797                              nest_offl%chem_top(0:1,nys:nyn,nxl:nxr,n),       &
798                              nxl+1, nys+1, nest_offl%tind+1,                  &
799                              nxr-nxl+1, nyn-nys+1, 2, .TRUE. )
800                nest_offl%chem_from_file_t(n) = .TRUE.
801             ENDIF
802          ENDDO
803       ENDIF
804
805!
806!--    Close input file
807       CALL close_input_file( pids_id )
808#endif
809!
810!--    Set control flag to indicate that boundary data has been initially
811!--    input.
812       nest_offl%init = .TRUE.
813!
814!--    Call offline nesting for salsa
815       IF ( salsa )  CALL salsa_nesting_offl_input
816!
817!--    End of CPU measurement
818       CALL cpu_log( log_point_s(86), 'NetCDF input forcing', 'stop' )
819
820    END SUBROUTINE nesting_offl_input
821
822
823!------------------------------------------------------------------------------!
824! Description:
825! ------------
826!> In this subroutine a constant mass within the model domain is guaranteed.
827!> Larger-scale models may be based on a compressible equation system, which is
828!> not consistent with PALMs incompressible equation system. In order to avoid
829!> a decrease or increase of mass during the simulation, non-divergent flow
830!> through the lateral and top boundaries is compensated by the vertical wind
831!> component at the top boundary.
832!------------------------------------------------------------------------------!
833    SUBROUTINE nesting_offl_mass_conservation
834
835       INTEGER(iwp) ::  i !< grid index in x-direction
836       INTEGER(iwp) ::  j !< grid index in y-direction
837       INTEGER(iwp) ::  k !< grid index in z-direction
838
839       REAL(wp) ::  d_area_t                        !< inverse of the total area of the horizontal model domain
840       REAL(wp) ::  w_correct                       !< vertical velocity increment required to compensate non-divergent flow through the boundaries
841       REAL(wp), DIMENSION(1:3) ::  volume_flow_l   !< local volume flow
842
843
844       IF ( debug_output_timestep )  CALL debug_message( 'nesting_offl_mass_conservation', 'start' )
845
846       CALL  cpu_log( log_point(58), 'offline nesting', 'start' )
847       
848       volume_flow   = 0.0_wp
849       volume_flow_l = 0.0_wp
850
851       d_area_t = 1.0_wp / ( ( nx + 1 ) * dx * ( ny + 1 ) * dy )
852
853       IF ( bc_dirichlet_l )  THEN
854          i = nxl
855          DO  j = nys, nyn
856             DO  k = nzb+1, nzt
857                volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k) * dy   &
858                                     * MERGE( 1.0_wp, 0.0_wp,                  &
859                                              BTEST( wall_flags_0(k,j,i), 1 ) )
860             ENDDO
861          ENDDO
862       ENDIF
863       IF ( bc_dirichlet_r )  THEN
864          i = nxr+1
865          DO  j = nys, nyn
866             DO  k = nzb+1, nzt
867                volume_flow_l(1) = volume_flow_l(1) - u(k,j,i) * dzw(k) * dy   &
868                                     * MERGE( 1.0_wp, 0.0_wp,                  &
869                                              BTEST( wall_flags_0(k,j,i), 1 ) )
870             ENDDO
871          ENDDO
872       ENDIF
873       IF ( bc_dirichlet_s )  THEN
874          j = nys
875          DO  i = nxl, nxr
876             DO  k = nzb+1, nzt
877                volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k) * dx   &
878                                     * MERGE( 1.0_wp, 0.0_wp,                  &
879                                              BTEST( wall_flags_0(k,j,i), 2 ) )
880             ENDDO
881          ENDDO
882       ENDIF
883       IF ( bc_dirichlet_n )  THEN
884          j = nyn+1
885          DO  i = nxl, nxr
886             DO  k = nzb+1, nzt
887                volume_flow_l(2) = volume_flow_l(2) - v(k,j,i) * dzw(k) * dx   &
888                                     * MERGE( 1.0_wp, 0.0_wp,                  &
889                                              BTEST( wall_flags_0(k,j,i), 2 ) )
890             ENDDO
891          ENDDO
892       ENDIF
893!
894!--    Top boundary
895       k = nzt
896       DO  i = nxl, nxr
897          DO  j = nys, nyn
898             volume_flow_l(3) = volume_flow_l(3) - w(k,j,i) * dx * dy
899          ENDDO
900       ENDDO
901
902#if defined( __parallel )
903       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
904       CALL MPI_ALLREDUCE( volume_flow_l, volume_flow, 3, MPI_REAL, MPI_SUM,   &
905                           comm2d, ierr )
906#else
907       volume_flow = volume_flow_l
908#endif
909
910       w_correct = SUM( volume_flow ) * d_area_t
911
912       DO  i = nxl, nxr
913          DO  j = nys, nyn
914             DO  k = nzt, nzt + 1
915                w(k,j,i) = w(k,j,i) + w_correct
916             ENDDO
917          ENDDO
918       ENDDO
919       
920       CALL  cpu_log( log_point(58), 'offline nesting', 'stop' )
921
922       IF ( debug_output_timestep )  CALL debug_message( 'nesting_offl_mass_conservation', 'end' )
923
924    END SUBROUTINE nesting_offl_mass_conservation
925
926
927!------------------------------------------------------------------------------!
928! Description:
929! ------------
930!> Set the lateral and top boundary conditions in case the PALM domain is
931!> nested offline in a mesoscale model. Further, average boundary data and
932!> determine mean profiles, further used for correct damping in the sponge
933!> layer.
934!------------------------------------------------------------------------------!
935    SUBROUTINE nesting_offl_bc                     
936
937       INTEGER(iwp) ::  i !< running index x-direction
938       INTEGER(iwp) ::  j !< running index y-direction
939       INTEGER(iwp) ::  k !< running index z-direction
940       INTEGER(iwp) ::  n !< running index for chemical species
941       
942       REAL(wp), DIMENSION(nzb:nzt+1) ::  pt_ref   !< reference profile for potential temperature
943       REAL(wp), DIMENSION(nzb:nzt+1) ::  pt_ref_l !< reference profile for potential temperature on subdomain
944       REAL(wp), DIMENSION(nzb:nzt+1) ::  q_ref    !< reference profile for mixing ratio
945       REAL(wp), DIMENSION(nzb:nzt+1) ::  q_ref_l  !< reference profile for mixing ratio on subdomain
946       REAL(wp), DIMENSION(nzb:nzt+1) ::  u_ref    !< reference profile for u-component
947       REAL(wp), DIMENSION(nzb:nzt+1) ::  u_ref_l  !< reference profile for u-component on subdomain
948       REAL(wp), DIMENSION(nzb:nzt+1) ::  v_ref    !< reference profile for v-component
949       REAL(wp), DIMENSION(nzb:nzt+1) ::  v_ref_l  !< reference profile for v-component on subdomain
950
951       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ref_chem   !< reference profile for chemical species
952       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ref_chem_l !< reference profile for chemical species on subdomain
953
954       IF ( debug_output_timestep )  CALL debug_message( 'nesting_offl_bc', 'start' )
955
956       CALL  cpu_log( log_point(58), 'offline nesting', 'start' )     
957!
958!--    Initialize mean profiles, derived from boundary data, to zero
959       pt_ref   = 0.0_wp
960       q_ref    = 0.0_wp
961       u_ref    = 0.0_wp
962       v_ref    = 0.0_wp
963       
964       pt_ref_l = 0.0_wp
965       q_ref_l  = 0.0_wp
966       u_ref_l  = 0.0_wp
967       v_ref_l  = 0.0_wp
968!
969!--    If required, allocate temporary arrays to compute chemistry mean profiles
970       IF ( air_chemistry )  THEN
971          ALLOCATE( ref_chem(nzb:nzt+1,1:UBOUND( chem_species, 1 ) )   )
972          ALLOCATE( ref_chem_l(nzb:nzt+1,1:UBOUND( chem_species, 1 ) ) )
973          ref_chem   = 0.0_wp
974          ref_chem_l = 0.0_wp
975       ENDIF
976!
977!--    Set boundary conditions of u-, v-, w-component, as well as q, and pt.
978!--    Note, boundary values at the left boundary: i=-1 (v,w,pt,q) and
979!--    i=0 (u), at the right boundary: i=nxr+1 (all), at the south boundary:
980!--    j=-1 (u,w,pt,q) and j=0 (v), at the north boundary: j=nyn+1 (all).
981!--    Please note, at the left (for u) and south (for v) boundary, values
982!--    for u and v are set also at i/j=-1, since these values are used in
983!--    boundary_conditions() to restore prognostic values.
984!--    Further, sum up data to calculate mean profiles from boundary data,
985!--    used for Rayleigh damping.
986       IF ( bc_dirichlet_l )  THEN
987
988          DO  j = nys, nyn
989             DO  k = nzb+1, nzt
990                u(k,j,0) = interpolate_in_time( nest_offl%u_left(0,k,j),       &
991                                                nest_offl%u_left(1,k,j),       &
992                                                fac_dt ) *                     &
993                             MERGE( 1.0_wp, 0.0_wp,                            &
994                                    BTEST( wall_flags_0(k,j,0), 1 ) )
995                u(k,j,-1) = u(k,j,0)
996             ENDDO
997             u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,j,0)
998          ENDDO
999
1000          DO  j = nys, nyn
1001             DO  k = nzb+1, nzt-1
1002                w(k,j,-1) = interpolate_in_time( nest_offl%w_left(0,k,j),      &
1003                                                 nest_offl%w_left(1,k,j),      &
1004                                                 fac_dt ) *                    &
1005                            MERGE( 1.0_wp, 0.0_wp,                             &
1006                                   BTEST( wall_flags_0(k,j,-1), 3 ) )
1007             ENDDO
1008             w(nzt,j,-1) = w(nzt-1,j,-1)
1009          ENDDO
1010
1011          DO  j = nysv, nyn
1012             DO  k = nzb+1, nzt
1013                v(k,j,-1) = interpolate_in_time( nest_offl%v_left(0,k,j),      &
1014                                                 nest_offl%v_left(1,k,j),      &
1015                                                 fac_dt ) *                    &
1016                               MERGE( 1.0_wp, 0.0_wp,                          &
1017                                      BTEST( wall_flags_0(k,j,-1), 2 ) )
1018             ENDDO
1019             v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,j,-1)
1020          ENDDO
1021
1022          IF ( .NOT. neutral )  THEN
1023             DO  j = nys, nyn
1024                DO  k = nzb+1, nzt
1025                   pt(k,j,-1) = interpolate_in_time( nest_offl%pt_left(0,k,j), &
1026                                                     nest_offl%pt_left(1,k,j), &
1027                                                     fac_dt )
1028 
1029                ENDDO
1030                pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,j,-1)
1031             ENDDO
1032          ENDIF
1033
1034          IF ( humidity )  THEN
1035             DO  j = nys, nyn
1036                DO  k = nzb+1, nzt
1037                   q(k,j,-1) = interpolate_in_time( nest_offl%q_left(0,k,j),   &
1038                                                    nest_offl%q_left(1,k,j),   &
1039                                                    fac_dt )
1040 
1041                ENDDO
1042                q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,j,-1)
1043             ENDDO
1044          ENDIF
1045         
1046          IF ( air_chemistry )  THEN
1047             DO  n = 1, UBOUND( chem_species, 1 )
1048                IF ( nest_offl%chem_from_file_l(n) )  THEN                   
1049                   DO  j = nys, nyn
1050                      DO  k = nzb+1, nzt
1051                         chem_species(n)%conc(k,j,-1) = interpolate_in_time(   &
1052                                                  nest_offl%chem_left(0,k,j,n),&
1053                                                  nest_offl%chem_left(1,k,j,n),&
1054                                                  fac_dt                   )
1055                      ENDDO
1056                      ref_chem_l(nzb+1:nzt,n) = ref_chem_l(nzb+1:nzt,n)        &
1057                                         + chem_species(n)%conc(nzb+1:nzt,j,-1)
1058                   ENDDO
1059                ENDIF
1060             ENDDO
1061          ENDIF
1062
1063       ENDIF
1064
1065       IF ( bc_dirichlet_r )  THEN
1066
1067          DO  j = nys, nyn
1068             DO  k = nzb+1, nzt
1069                u(k,j,nxr+1) = interpolate_in_time( nest_offl%u_right(0,k,j),  &
1070                                                    nest_offl%u_right(1,k,j),  &
1071                                                    fac_dt ) *                 &
1072                             MERGE( 1.0_wp, 0.0_wp,                            &
1073                                    BTEST( wall_flags_0(k,j,nxr+1), 1 ) )
1074             ENDDO
1075             u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,j,nxr+1)
1076          ENDDO
1077          DO  j = nys, nyn
1078             DO  k = nzb+1, nzt-1
1079                w(k,j,nxr+1) = interpolate_in_time( nest_offl%w_right(0,k,j),  &
1080                                                    nest_offl%w_right(1,k,j),  &
1081                                                    fac_dt ) *                 &
1082                             MERGE( 1.0_wp, 0.0_wp,                            &
1083                                    BTEST( wall_flags_0(k,j,nxr+1), 3 ) )
1084             ENDDO
1085             w(nzt,j,nxr+1) = w(nzt-1,j,nxr+1)
1086          ENDDO
1087
1088          DO  j = nysv, nyn
1089             DO  k = nzb+1, nzt
1090                v(k,j,nxr+1) = interpolate_in_time( nest_offl%v_right(0,k,j),  &
1091                                                    nest_offl%v_right(1,k,j),  &
1092                                                    fac_dt ) *                 &
1093                             MERGE( 1.0_wp, 0.0_wp,                            &
1094                                    BTEST( wall_flags_0(k,j,nxr+1), 2 ) )
1095             ENDDO
1096             v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,j,nxr+1)
1097          ENDDO
1098
1099          IF ( .NOT. neutral )  THEN
1100             DO  j = nys, nyn
1101                DO  k = nzb+1, nzt
1102                   pt(k,j,nxr+1) = interpolate_in_time(                        &
1103                                                  nest_offl%pt_right(0,k,j),   &
1104                                                  nest_offl%pt_right(1,k,j),   &
1105                                                  fac_dt )
1106                ENDDO
1107                pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,j,nxr+1)
1108             ENDDO
1109          ENDIF
1110
1111          IF ( humidity )  THEN
1112             DO  j = nys, nyn
1113                DO  k = nzb+1, nzt
1114                   q(k,j,nxr+1) = interpolate_in_time(                         &
1115                                                  nest_offl%q_right(0,k,j),    &
1116                                                  nest_offl%q_right(1,k,j),    &
1117                                                  fac_dt ) 
1118 
1119                ENDDO
1120                q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,j,nxr+1)
1121             ENDDO
1122          ENDIF
1123         
1124          IF ( air_chemistry )  THEN
1125             DO  n = 1, UBOUND( chem_species, 1 )
1126                IF ( nest_offl%chem_from_file_r(n) )  THEN 
1127                   DO  j = nys, nyn
1128                      DO  k = nzb+1, nzt
1129                         chem_species(n)%conc(k,j,nxr+1) = interpolate_in_time(&
1130                                                 nest_offl%chem_right(0,k,j,n),&
1131                                                 nest_offl%chem_right(1,k,j,n),&
1132                                                 fac_dt                       )
1133                      ENDDO
1134                      ref_chem_l(nzb+1:nzt,n) = ref_chem_l(nzb+1:nzt,n)        &
1135                                       + chem_species(n)%conc(nzb+1:nzt,j,nxr+1)
1136                   ENDDO
1137                ENDIF
1138             ENDDO
1139          ENDIF
1140
1141       ENDIF
1142
1143       IF ( bc_dirichlet_s )  THEN
1144
1145          DO  i = nxl, nxr
1146             DO  k = nzb+1, nzt
1147                v(k,0,i) = interpolate_in_time( nest_offl%v_south(0,k,i),      &
1148                                                nest_offl%v_south(1,k,i),      &
1149                                                fac_dt ) *                     &
1150                           MERGE( 1.0_wp, 0.0_wp,                              &
1151                                  BTEST( wall_flags_0(k,0,i), 2 ) )
1152                v(k,-1,i) = v(k,0,i) 
1153             ENDDO
1154             v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,0,i)
1155          ENDDO
1156
1157          DO  i = nxl, nxr
1158             DO  k = nzb+1, nzt-1
1159                w(k,-1,i) = interpolate_in_time( nest_offl%w_south(0,k,i),     &
1160                                                 nest_offl%w_south(1,k,i),     &
1161                                                 fac_dt ) *                    &
1162                           MERGE( 1.0_wp, 0.0_wp,                              &
1163                                  BTEST( wall_flags_0(k,-1,i), 3 ) )
1164             ENDDO
1165             w(nzt,-1,i) = w(nzt-1,-1,i)
1166          ENDDO
1167
1168          DO  i = nxlu, nxr
1169             DO  k = nzb+1, nzt
1170                u(k,-1,i) = interpolate_in_time( nest_offl%u_south(0,k,i),     &
1171                                                 nest_offl%u_south(1,k,i),     &
1172                                                 fac_dt ) *                    &
1173                           MERGE( 1.0_wp, 0.0_wp,                              &
1174                                  BTEST( wall_flags_0(k,-1,i), 1 ) )
1175             ENDDO
1176             u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,-1,i)
1177          ENDDO
1178
1179          IF ( .NOT. neutral )  THEN
1180             DO  i = nxl, nxr
1181                DO  k = nzb+1, nzt
1182                   pt(k,-1,i) = interpolate_in_time(                           &
1183                                                 nest_offl%pt_south(0,k,i),    &
1184                                                 nest_offl%pt_south(1,k,i),    &
1185                                                 fac_dt )
1186 
1187                ENDDO
1188                pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,-1,i)
1189             ENDDO
1190          ENDIF
1191
1192          IF ( humidity )  THEN
1193             DO  i = nxl, nxr
1194                DO  k = nzb+1, nzt
1195                   q(k,-1,i) = interpolate_in_time(                            &
1196                                                 nest_offl%q_south(0,k,i),     &
1197                                                 nest_offl%q_south(1,k,i),     &
1198                                                 fac_dt )
1199 
1200                ENDDO
1201                q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,-1,i)
1202             ENDDO
1203          ENDIF
1204         
1205          IF ( air_chemistry )  THEN
1206             DO  n = 1, UBOUND( chem_species, 1 )
1207                IF ( nest_offl%chem_from_file_s(n) )  THEN 
1208                   DO  i = nxl, nxr
1209                      DO  k = nzb+1, nzt
1210                         chem_species(n)%conc(k,-1,i) = interpolate_in_time(   &
1211                                                 nest_offl%chem_south(0,k,i,n),&
1212                                                 nest_offl%chem_south(1,k,i,n),&
1213                                                 fac_dt                    )
1214                      ENDDO
1215                      ref_chem_l(nzb+1:nzt,n) = ref_chem_l(nzb+1:nzt,n)        &
1216                                       + chem_species(n)%conc(nzb+1:nzt,-1,i)
1217                   ENDDO
1218                ENDIF
1219             ENDDO
1220          ENDIF
1221
1222       ENDIF
1223
1224       IF ( bc_dirichlet_n )  THEN
1225
1226          DO  i = nxl, nxr
1227             DO  k = nzb+1, nzt
1228                v(k,nyn+1,i) = interpolate_in_time( nest_offl%v_north(0,k,i),  &
1229                                                    nest_offl%v_north(1,k,i),  &
1230                                                    fac_dt ) *                 &
1231                               MERGE( 1.0_wp, 0.0_wp,                          &
1232                                      BTEST( wall_flags_0(k,nyn+1,i), 2 ) )
1233             ENDDO
1234             v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,nyn+1,i)
1235          ENDDO
1236          DO  i = nxl, nxr
1237             DO  k = nzb+1, nzt-1
1238                w(k,nyn+1,i) = interpolate_in_time( nest_offl%w_north(0,k,i),  &
1239                                                    nest_offl%w_north(1,k,i),  &
1240                                                    fac_dt ) *                 &
1241                               MERGE( 1.0_wp, 0.0_wp,                          &
1242                                      BTEST( wall_flags_0(k,nyn+1,i), 3 ) )
1243             ENDDO
1244             w(nzt,nyn+1,i) = w(nzt-1,nyn+1,i)
1245          ENDDO
1246
1247          DO  i = nxlu, nxr
1248             DO  k = nzb+1, nzt
1249                u(k,nyn+1,i) = interpolate_in_time( nest_offl%u_north(0,k,i),  &
1250                                                    nest_offl%u_north(1,k,i),  &
1251                                                    fac_dt ) *                 &
1252                               MERGE( 1.0_wp, 0.0_wp,                          &
1253                                      BTEST( wall_flags_0(k,nyn+1,i), 1 ) )
1254
1255             ENDDO
1256             u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,nyn+1,i)
1257          ENDDO
1258
1259          IF ( .NOT. neutral )  THEN
1260             DO  i = nxl, nxr
1261                DO  k = nzb+1, nzt
1262                   pt(k,nyn+1,i) = interpolate_in_time(                        &
1263                                                    nest_offl%pt_north(0,k,i), &
1264                                                    nest_offl%pt_north(1,k,i), &
1265                                                    fac_dt )
1266 
1267                ENDDO
1268                pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,nyn+1,i)
1269             ENDDO
1270          ENDIF
1271
1272          IF ( humidity )  THEN
1273             DO  i = nxl, nxr
1274                DO  k = nzb+1, nzt
1275                   q(k,nyn+1,i) = interpolate_in_time(                         &
1276                                                    nest_offl%q_north(0,k,i),  &
1277                                                    nest_offl%q_north(1,k,i),  &
1278                                                    fac_dt )
1279 
1280                ENDDO
1281                q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,nyn+1,i)
1282             ENDDO
1283          ENDIF
1284         
1285          IF ( air_chemistry )  THEN
1286             DO  n = 1, UBOUND( chem_species, 1 )
1287                IF ( nest_offl%chem_from_file_n(n) )  THEN 
1288                   DO  i = nxl, nxr
1289                      DO  k = nzb+1, nzt
1290                         chem_species(n)%conc(k,nyn+1,i) = interpolate_in_time(&
1291                                                 nest_offl%chem_north(0,k,i,n),&
1292                                                 nest_offl%chem_north(1,k,i,n),&
1293                                                 fac_dt                       )
1294                      ENDDO
1295                      ref_chem_l(nzb+1:nzt,n) = ref_chem_l(nzb+1:nzt,n)        &
1296                                       + chem_species(n)%conc(nzb+1:nzt,nyn+1,i)
1297                   ENDDO
1298                ENDIF
1299             ENDDO
1300          ENDIF
1301
1302       ENDIF
1303!
1304!--    Top boundary
1305       DO  i = nxlu, nxr
1306          DO  j = nys, nyn
1307             u(nzt+1,j,i) = interpolate_in_time( nest_offl%u_top(0,j,i),       &
1308                                                 nest_offl%u_top(1,j,i),       &
1309                                                 fac_dt ) *                    &
1310                           MERGE( 1.0_wp, 0.0_wp,                              &
1311                                  BTEST( wall_flags_0(nzt+1,j,i), 1 ) )
1312             u_ref_l(nzt+1) = u_ref_l(nzt+1) + u(nzt+1,j,i)
1313          ENDDO
1314       ENDDO
1315!
1316!--    For left boundary set boundary condition for u-component also at top
1317!--    grid point.
1318!--    Note, this has no effect on the numeric solution, only for data output.
1319       IF ( bc_dirichlet_l )  u(nzt+1,:,nxl) = u(nzt+1,:,nxlu)
1320
1321       DO  i = nxl, nxr
1322          DO  j = nysv, nyn
1323             v(nzt+1,j,i) = interpolate_in_time( nest_offl%v_top(0,j,i),       &
1324                                                 nest_offl%v_top(1,j,i),       &
1325                                                 fac_dt ) *                    &
1326                           MERGE( 1.0_wp, 0.0_wp,                              &
1327                                  BTEST( wall_flags_0(nzt+1,j,i), 2 ) )
1328             v_ref_l(nzt+1) = v_ref_l(nzt+1) + v(nzt+1,j,i)
1329          ENDDO
1330       ENDDO
1331!
1332!--    For south boundary set boundary condition for v-component also at top
1333!--    grid point.
1334!--    Note, this has no effect on the numeric solution, only for data output.
1335       IF ( bc_dirichlet_s )  v(nzt+1,nys,:) = v(nzt+1,nysv,:)
1336
1337       DO  i = nxl, nxr
1338          DO  j = nys, nyn
1339             w(nzt,j,i) = interpolate_in_time( nest_offl%w_top(0,j,i),         &
1340                                               nest_offl%w_top(1,j,i),         &
1341                                               fac_dt ) *                      &
1342                           MERGE( 1.0_wp, 0.0_wp,                              &
1343                                  BTEST( wall_flags_0(nzt,j,i), 3 ) )
1344             w(nzt+1,j,i) = w(nzt,j,i)
1345          ENDDO
1346       ENDDO
1347
1348
1349       IF ( .NOT. neutral )  THEN
1350          DO  i = nxl, nxr
1351             DO  j = nys, nyn
1352                pt(nzt+1,j,i) = interpolate_in_time( nest_offl%pt_top(0,j,i),  &
1353                                                     nest_offl%pt_top(1,j,i),  &
1354                                                     fac_dt )
1355                pt_ref_l(nzt+1) = pt_ref_l(nzt+1) + pt(nzt+1,j,i)
1356             ENDDO
1357          ENDDO
1358       ENDIF
1359
1360       IF ( humidity )  THEN
1361          DO  i = nxl, nxr
1362             DO  j = nys, nyn
1363                q(nzt+1,j,i) = interpolate_in_time( nest_offl%q_top(0,j,i),    &
1364                                                    nest_offl%q_top(1,j,i),    &
1365                                                    fac_dt )
1366                q_ref_l(nzt+1) = q_ref_l(nzt+1) + q(nzt+1,j,i)
1367             ENDDO
1368          ENDDO
1369       ENDIF
1370       
1371       IF ( air_chemistry )  THEN
1372          DO  n = 1, UBOUND( chem_species, 1 )
1373             IF ( nest_offl%chem_from_file_t(n) )  THEN 
1374                DO  i = nxl, nxr
1375                   DO  j = nys, nyn
1376                      chem_species(n)%conc(nzt+1,j,i) = interpolate_in_time(   &
1377                                              nest_offl%chem_top(0,j,i,n),     &
1378                                              nest_offl%chem_top(1,j,i,n),     &
1379                                              fac_dt                       )
1380                      ref_chem_l(nzt+1,n) = ref_chem_l(nzt+1,n) +              &
1381                                            chem_species(n)%conc(nzt+1,j,i)
1382                   ENDDO
1383                ENDDO
1384             ENDIF
1385          ENDDO
1386       ENDIF
1387!
1388!--    Moreover, set Neumann boundary condition for subgrid-scale TKE,
1389!--    passive scalar, dissipation, and chemical species if required
1390       IF ( rans_mode  .AND.  rans_tke_e )  THEN
1391          IF (  bc_dirichlet_l )  diss(:,:,nxl-1) = diss(:,:,nxl)
1392          IF (  bc_dirichlet_r )  diss(:,:,nxr+1) = diss(:,:,nxr)
1393          IF (  bc_dirichlet_s )  diss(:,nys-1,:) = diss(:,nys,:)
1394          IF (  bc_dirichlet_n )  diss(:,nyn+1,:) = diss(:,nyn,:)
1395       ENDIF
1396!        IF ( .NOT. constant_diffusion )  THEN
1397!           IF (  bc_dirichlet_l )  e(:,:,nxl-1) = e(:,:,nxl)
1398!           IF (  bc_dirichlet_r )  e(:,:,nxr+1) = e(:,:,nxr)
1399!           IF (  bc_dirichlet_s )  e(:,nys-1,:) = e(:,nys,:)
1400!           IF (  bc_dirichlet_n )  e(:,nyn+1,:) = e(:,nyn,:)
1401!           e(nzt+1,:,:) = e(nzt,:,:)
1402!        ENDIF
1403!        IF ( passive_scalar )  THEN
1404!           IF (  bc_dirichlet_l )  s(:,:,nxl-1) = s(:,:,nxl)
1405!           IF (  bc_dirichlet_r )  s(:,:,nxr+1) = s(:,:,nxr)
1406!           IF (  bc_dirichlet_s )  s(:,nys-1,:) = s(:,nys,:)
1407!           IF (  bc_dirichlet_n )  s(:,nyn+1,:) = s(:,nyn,:)
1408!        ENDIF
1409
1410       CALL exchange_horiz( u, nbgp )
1411       CALL exchange_horiz( v, nbgp )
1412       CALL exchange_horiz( w, nbgp )
1413       IF ( .NOT. neutral )  CALL exchange_horiz( pt, nbgp )
1414       IF ( humidity      )  CALL exchange_horiz( q,  nbgp )
1415       IF ( air_chemistry )  THEN
1416          DO  n = 1, UBOUND( chem_species, 1 )
1417!
1418!--          Do local exchange only when necessary, i.e. when data is coming
1419!--          from dynamic file.
1420             IF ( nest_offl%chem_from_file_t(n) )                              &
1421                CALL exchange_horiz( chem_species(n)%conc, nbgp )
1422          ENDDO
1423       ENDIF
1424!
1425!--    Set top boundary condition at all horizontal grid points, also at the
1426!--    lateral boundary grid points.
1427       w(nzt+1,:,:) = w(nzt,:,:)       
1428!
1429!--    Offline nesting for salsa
1430       IF ( salsa )  CALL salsa_nesting_offl_bc
1431!
1432!--    In case of Rayleigh damping, where the profiles u_init, v_init
1433!--    q_init and pt_init are still used, update these profiles from the
1434!--    averaged boundary data.
1435!--    But first, average these data.
1436#if defined( __parallel )
1437       CALL MPI_ALLREDUCE( u_ref_l, u_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM,     &
1438                           comm2d, ierr )
1439       CALL MPI_ALLREDUCE( v_ref_l, v_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM,     &
1440                           comm2d, ierr )
1441       IF ( humidity )  THEN
1442          CALL MPI_ALLREDUCE( q_ref_l, q_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM,  &
1443                              comm2d, ierr )
1444       ENDIF
1445       IF ( .NOT. neutral )  THEN
1446          CALL MPI_ALLREDUCE( pt_ref_l, pt_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM,&
1447                              comm2d, ierr )
1448       ENDIF
1449       IF ( air_chemistry )  THEN
1450          CALL MPI_ALLREDUCE( ref_chem_l, ref_chem,                            &
1451                              ( nzt+1-nzb+1 ) * SIZE( ref_chem(nzb,:) ),       &
1452                              MPI_REAL, MPI_SUM, comm2d, ierr )
1453       ENDIF
1454#else
1455       u_ref  = u_ref_l
1456       v_ref  = v_ref_l
1457       IF ( humidity )       q_ref    = q_ref_l
1458       IF ( .NOT. neutral )  pt_ref   = pt_ref_l
1459       IF ( air_chemistry )  ref_chem = ref_chem_l
1460#endif
1461!
1462!--    Average data. Note, reference profiles up to nzt are derived from lateral
1463!--    boundaries, at the model top it is derived from the top boundary. Thus,
1464!--    number of input data is different from nzb:nzt compared to nzt+1.
1465!--    Derived from lateral boundaries.
1466       u_ref(nzb:nzt) = u_ref(nzb:nzt) / REAL( 2.0_wp * ( ny + 1 + nx     ),   &
1467                                               KIND = wp ) 
1468       v_ref(nzb:nzt) = v_ref(nzb:nzt) / REAL( 2.0_wp * ( ny   + nx + 1   ),   &
1469                                               KIND = wp )
1470       IF ( humidity )                                                         &
1471          q_ref(nzb:nzt) = q_ref(nzb:nzt)   / REAL( 2.0_wp *                   &
1472                                                          ( ny + 1 + nx + 1 ), &
1473                                                    KIND = wp )
1474       IF ( .NOT. neutral )                                                    &
1475          pt_ref(nzb:nzt) = pt_ref(nzb:nzt) / REAL( 2.0_wp *                   &
1476                                                          ( ny + 1 + nx + 1 ), &
1477                                              KIND = wp )
1478       IF ( air_chemistry )                                                    &
1479          ref_chem(nzb:nzt,:) = ref_chem(nzb:nzt,:) / REAL( 2.0_wp *           &
1480                                                          ( ny + 1 + nx + 1 ), &
1481                                                            KIND = wp )
1482!
1483!--    Derived from top boundary.   
1484       u_ref(nzt+1) = u_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx     ), KIND = wp ) 
1485       v_ref(nzt+1) = v_ref(nzt+1) / REAL( ( ny     ) * ( nx + 1 ), KIND = wp )
1486       IF ( humidity )                                                         &
1487          q_ref(nzt+1) = q_ref(nzt+1)   / REAL( ( ny + 1 ) * ( nx + 1 ),       &
1488                                                KIND = wp )
1489       IF ( .NOT. neutral )                                                    &
1490          pt_ref(nzt+1) = pt_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx + 1 ),       &
1491                                                KIND = wp )
1492       IF ( air_chemistry )                                                    &
1493          ref_chem(nzt+1,:) = ref_chem(nzt+1,:) /                              &
1494                              REAL( ( ny + 1 ) * ( nx + 1 ),KIND = wp )
1495!
1496!--    Write onto init profiles, which are used for damping. Also set lower
1497!--    boundary condition for scalars (not required for u and v as these are
1498!--    zero at k=nzb.
1499       u_init = u_ref
1500       v_init = v_ref
1501       IF ( humidity      )  THEN
1502          q_init      = q_ref
1503          q_init(nzb) = q_init(nzb+1)
1504       ENDIF
1505       IF ( .NOT. neutral )  THEN
1506          pt_init      = pt_ref
1507          pt_init(nzb) = pt_init(nzb+1)
1508       ENDIF
1509
1510       IF ( air_chemistry )  THEN
1511          DO  n = 1, UBOUND( chem_species, 1 )
1512             IF ( nest_offl%chem_from_file_t(n) )  THEN
1513                chem_species(n)%conc_pr_init(:) = ref_chem(:,n)
1514                chem_species(n)%conc_pr_init(nzb) =                            &
1515                                            chem_species(n)%conc_pr_init(nzb+1)
1516             ENDIF
1517          ENDDO
1518       ENDIF
1519
1520       IF ( ALLOCATED( ref_chem   ) )  DEALLOCATE( ref_chem   )
1521       IF ( ALLOCATED( ref_chem_l ) )  DEALLOCATE( ref_chem_l )     
1522!
1523!--    Further, adjust Rayleigh damping height in case of time-changing conditions.
1524!--    Therefore, calculate boundary-layer depth first.
1525       CALL nesting_offl_calc_zi
1526       CALL adjust_sponge_layer 
1527       
1528       CALL  cpu_log( log_point(58), 'offline nesting', 'stop' )
1529
1530       IF ( debug_output_timestep )  CALL debug_message( 'nesting_offl_bc', 'end' )
1531
1532
1533    END SUBROUTINE nesting_offl_bc
1534
1535!------------------------------------------------------------------------------!
1536! Description:
1537!------------------------------------------------------------------------------!
1538!>  Update of the geostrophic wind components.
1539!>  @todo: update geostrophic wind also in the child domains (should be done
1540!>         in the nesting.
1541!------------------------------------------------------------------------------!
1542    SUBROUTINE nesting_offl_geostrophic_wind
1543
1544       INTEGER(iwp) ::  k
1545!
1546!--    Update geostrophic wind components from dynamic input file.
1547       DO  k = nzb+1, nzt
1548          ug(k) = interpolate_in_time( nest_offl%ug(0,k), nest_offl%ug(1,k),   &
1549                                       fac_dt )
1550          vg(k) = interpolate_in_time( nest_offl%vg(0,k), nest_offl%vg(1,k),   &
1551                                       fac_dt )
1552       ENDDO
1553       ug(nzt+1) = ug(nzt)
1554       vg(nzt+1) = vg(nzt)
1555
1556    END SUBROUTINE nesting_offl_geostrophic_wind
1557
1558!------------------------------------------------------------------------------!
1559! Description:
1560!------------------------------------------------------------------------------!
1561!>  Determine the interpolation constant for time interpolation. The
1562!>  calculation is separated from the nesting_offl_bc and
1563!>  nesting_offl_geostrophic_wind in order to be independent on the order
1564!>  of calls.
1565!------------------------------------------------------------------------------!
1566    SUBROUTINE nesting_offl_interpolation_factor
1567!
1568!--    Determine interpolation factor and limit it to 1. This is because
1569!--    t+dt can slightly exceed time(tind_p) before boundary data is updated
1570!--    again.
1571       fac_dt = ( time_utc_init + time_since_reference_point                   &
1572                - nest_offl%time(nest_offl%tind) + dt_3d ) /                   &
1573           ( nest_offl%time(nest_offl%tind_p) - nest_offl%time(nest_offl%tind) )
1574
1575       fac_dt = MIN( 1.0_wp, fac_dt )
1576
1577    END SUBROUTINE nesting_offl_interpolation_factor
1578
1579!------------------------------------------------------------------------------!
1580! Description:
1581!------------------------------------------------------------------------------!
1582!> Calculates the boundary-layer depth from the boundary data, according to
1583!> bulk-Richardson criterion.
1584!------------------------------------------------------------------------------!
1585    SUBROUTINE nesting_offl_calc_zi
1586
1587       INTEGER(iwp) :: i                            !< loop index in x-direction
1588       INTEGER(iwp) :: j                            !< loop index in y-direction
1589       INTEGER(iwp) :: k                            !< loop index in z-direction
1590       INTEGER(iwp) :: k_max_loc                    !< index of maximum wind speed along z-direction
1591       INTEGER(iwp) :: k_surface                    !< topography top index in z-direction
1592       INTEGER(iwp) :: num_boundary_gp_non_cyclic   !< number of non-cyclic boundaries, used for averaging ABL depth
1593       INTEGER(iwp) :: num_boundary_gp_non_cyclic_l !< number of non-cyclic boundaries, used for averaging ABL depth
1594   
1595       REAL(wp) ::  ri_bulk                 !< bulk Richardson number
1596       REAL(wp) ::  ri_bulk_crit = 0.25_wp  !< critical bulk Richardson number
1597       REAL(wp) ::  topo_max                !< maximum topography level in model domain
1598       REAL(wp) ::  topo_max_l              !< maximum topography level in subdomain
1599       REAL(wp) ::  vpt_surface             !< near-surface virtual potential temperature
1600       REAL(wp) ::  zi_l                    !< mean boundary-layer depth on subdomain   
1601       REAL(wp) ::  zi_local                !< local boundary-layer depth 
1602   
1603       REAL(wp), DIMENSION(nzb:nzt+1) ::  vpt_col !< vertical profile of virtual potential temperature at (j,i)-grid point
1604       REAL(wp), DIMENSION(nzb:nzt+1) ::  uv_abs  !< vertical profile of horizontal wind speed at (j,i)-grid point
1605
1606       
1607!
1608!--    Calculate mean boundary-layer height from boundary data.
1609!--    Start with the left and right boundaries.
1610       zi_l      = 0.0_wp
1611       num_boundary_gp_non_cyclic_l = 0
1612       IF ( bc_dirichlet_l  .OR.  bc_dirichlet_r )  THEN
1613!
1614!--       Sum-up and store number of boundary grid points used for averaging
1615!--       ABL depth
1616          num_boundary_gp_non_cyclic_l = num_boundary_gp_non_cyclic_l +        &
1617                                         nxr - nxl + 1
1618!
1619!--       Determine index along x. Please note, index indicates boundary
1620!--       grid point for scalars.
1621          i = MERGE( -1, nxr + 1, bc_dirichlet_l )
1622   
1623          DO  j = nys, nyn
1624!
1625!--          Determine topography top index at current (j,i) index
1626             k_surface = topo_top_ind(j,i,0)
1627!
1628!--          Pre-compute surface virtual temperature. Therefore, use 2nd
1629!--          prognostic level according to Heinze et al. (2017).
1630             IF ( humidity )  THEN
1631                vpt_surface = pt(k_surface+2,j,i) *                            &
1632                            ( 1.0_wp + 0.61_wp * q(k_surface+2,j,i) )
1633                vpt_col     = pt(:,j,i) * ( 1.0_wp + 0.61_wp * q(:,j,i) )
1634             ELSE
1635                vpt_surface = pt(k_surface+2,j,i)
1636                vpt_col     = pt(:,j,i)
1637             ENDIF
1638!
1639!--          Calculate local boundary layer height from bulk Richardson number,
1640!--          i.e. the height where the bulk Richardson number exceeds its
1641!--          critical value of 0.25 (according to Heinze et al., 2017).
1642!--          Note, no interpolation of u- and v-component is made, as both
1643!--          are mainly mean inflow profiles with very small spatial variation.
1644!--          Add a safety factor in case the velocity term becomes zero. This
1645!--          may happen if overhanging 3D structures are directly located at
1646!--          the boundary, where velocity inside the building is zero
1647!--          (k_surface is the index of the lowest upward-facing surface).
1648             uv_abs(:) = SQRT( MERGE( u(:,j,i+1), u(:,j,i),                   &
1649                                      bc_dirichlet_l )**2 +                   &
1650                               v(:,j,i)**2 )
1651!
1652!--          Determine index of the maximum wind speed                               
1653             k_max_loc = MAXLOC( uv_abs(:), DIM = 1 ) - 1
1654
1655             zi_local = 0.0_wp
1656             DO  k = k_surface+1, nzt
1657                ri_bulk = zu(k) * g / vpt_surface *                            &
1658                          ( vpt_col(k) - vpt_surface ) /                       &
1659                          ( uv_abs(k) + 1E-5_wp ) 
1660!
1661!--             Check if critical Richardson number is exceeded. Further, check
1662!--             if there is a maxium in the wind profile in order to detect also
1663!--             ABL heights in the stable boundary layer.
1664                IF ( zi_local == 0.0_wp  .AND.                                 &
1665                     ( ri_bulk > ri_bulk_crit  .OR.  k == k_max_loc ) )        &
1666                   zi_local = zu(k)           
1667             ENDDO
1668!
1669!--          Assure that the minimum local boundary-layer depth is at least at
1670!--          the second vertical grid level.
1671             zi_l = zi_l + MAX( zi_local, zu(k_surface+2) )   
1672             
1673          ENDDO
1674       
1675       ENDIF
1676!
1677!--    Do the same at the north and south boundaries.
1678       IF ( bc_dirichlet_s  .OR.  bc_dirichlet_n )  THEN
1679       
1680          num_boundary_gp_non_cyclic_l = num_boundary_gp_non_cyclic_l +        &
1681                                         nxr - nxl + 1
1682       
1683          j = MERGE( -1, nyn + 1, bc_dirichlet_s )
1684       
1685          DO  i = nxl, nxr
1686             k_surface = topo_top_ind(j,i,0)
1687 
1688             IF ( humidity )  THEN
1689                vpt_surface = pt(k_surface+2,j,i) *                            &
1690                            ( 1.0_wp + 0.61_wp * q(k_surface+2,j,i) )
1691                vpt_col     = pt(:,j,i) * ( 1.0_wp + 0.61_wp * q(:,j,i) )
1692             ELSE
1693                vpt_surface = pt(k_surface+2,j,i)
1694                vpt_col  = pt(:,j,i)
1695             ENDIF
1696         
1697             uv_abs(:) = SQRT( u(:,j,i)**2 +                                   &
1698                               MERGE( v(:,j+1,i), v(:,j,i),                    &
1699                               bc_dirichlet_s )**2 )
1700!
1701!--          Determine index of the maximum wind speed
1702             k_max_loc = MAXLOC( uv_abs(:), DIM = 1 ) - 1
1703         
1704             zi_local = 0.0_wp
1705             DO  k = k_surface+1, nzt               
1706                ri_bulk = zu(k) * g / vpt_surface *                            &
1707                       ( vpt_col(k) - vpt_surface ) /                          &
1708                       ( uv_abs(k) + 1E-5_wp ) 
1709!
1710!--             Check if critical Richardson number is exceeded. Further, check
1711!--             if there is a maxium in the wind profile in order to detect also
1712!--             ABL heights in the stable boundary layer.                       
1713                IF ( zi_local == 0.0_wp  .AND.                                 &
1714                     ( ri_bulk > ri_bulk_crit  .OR.  k == k_max_loc ) )        &
1715                   zi_local = zu(k)   
1716             ENDDO
1717             zi_l = zi_l + MAX( zi_local, zu(k_surface+2) )   
1718         
1719          ENDDO
1720         
1721       ENDIF
1722   
1723#if defined( __parallel )
1724       CALL MPI_ALLREDUCE( zi_l, zi_ribulk, 1, MPI_REAL, MPI_SUM,              &
1725                           comm2d, ierr )
1726       CALL MPI_ALLREDUCE( num_boundary_gp_non_cyclic_l,                       &
1727                           num_boundary_gp_non_cyclic,                         &
1728                           1, MPI_INTEGER, MPI_SUM, comm2d, ierr )
1729#else
1730       zi_ribulk = zi_l
1731       num_boundary_gp_non_cyclic = num_boundary_gp_non_cyclic_l
1732#endif
1733       zi_ribulk = zi_ribulk / REAL( num_boundary_gp_non_cyclic, KIND = wp )
1734!
1735!--    Finally, check if boundary layer depth is not below the any topography.
1736!--    zi_ribulk will be used to adjust rayleigh damping height, i.e. the
1737!--    lower level of the sponge layer, as well as to adjust the synthetic
1738!--    turbulence generator accordingly. If Rayleigh damping would be applied
1739!--    near buildings, etc., this would spoil the simulation results.
1740       topo_max_l = zw(MAXVAL( topo_top_ind(nys:nyn,nxl:nxr,0) ))
1741       
1742#if defined( __parallel )
1743       CALL MPI_ALLREDUCE( topo_max_l, topo_max, 1, MPI_REAL, MPI_MAX,         &
1744                           comm2d, ierr )
1745#else
1746       topo_max     = topo_max_l
1747#endif
1748!        zi_ribulk = MAX( zi_ribulk, topo_max )
1749       
1750    END SUBROUTINE nesting_offl_calc_zi
1751   
1752   
1753!------------------------------------------------------------------------------!
1754! Description:
1755!------------------------------------------------------------------------------!
1756!> Adjust the height where the rayleigh damping starts, i.e. the lower level
1757!> of the sponge layer.
1758!------------------------------------------------------------------------------!
1759    SUBROUTINE adjust_sponge_layer
1760
1761       INTEGER(iwp) :: k   !< loop index in z-direction
1762   
1763       REAL(wp) ::  rdh    !< updated Rayleigh damping height
1764 
1765   
1766       IF ( rayleigh_damping_height > 0.0_wp  .AND.                            &
1767            rayleigh_damping_factor > 0.0_wp )  THEN
1768!
1769!--       Update Rayleigh-damping height and re-calculate height-depending
1770!--       damping coefficients.
1771!--       Assure that rayleigh damping starts well above the boundary layer.   
1772          rdh = MIN( MAX( zi_ribulk * 1.3_wp, 10.0_wp * dz(1) ),               & 
1773                     0.8_wp * zu(nzt), rayleigh_damping_height )
1774!
1775!--       Update Rayleigh damping factor
1776          DO  k = nzb+1, nzt
1777             IF ( zu(k) >= rdh )  THEN
1778                rdf(k) = rayleigh_damping_factor *                             &
1779                                          ( SIN( pi * 0.5_wp * ( zu(k) - rdh ) &
1780                                        / ( zu(nzt) - rdh ) )                  &
1781                                          )**2
1782             ENDIF
1783          ENDDO
1784          rdf_sc = rdf
1785       
1786       ENDIF
1787
1788    END SUBROUTINE adjust_sponge_layer
1789   
1790!------------------------------------------------------------------------------!
1791! Description:
1792! ------------
1793!> Performs consistency checks
1794!------------------------------------------------------------------------------!
1795    SUBROUTINE nesting_offl_check_parameters 
1796!
1797!--    Check if offline nesting is applied in nested child domain.
1798       IF ( nesting_offline  .AND.  child_domain )  THEN
1799          message_string = 'Offline nesting is only applicable in root model.'
1800          CALL message( 'offline_nesting_check_parameters', 'PA0622', 1, 2, 0, 6, 0 )       
1801       ENDIF
1802
1803    END SUBROUTINE nesting_offl_check_parameters
1804   
1805!------------------------------------------------------------------------------!
1806! Description:
1807! ------------
1808!> Reads the parameter list nesting_offl_parameters
1809!------------------------------------------------------------------------------!
1810    SUBROUTINE nesting_offl_parin 
1811       
1812       CHARACTER (LEN=80) ::  line   !< dummy string that contains the current line of the parameter file
1813
1814
1815       NAMELIST /nesting_offl_parameters/   nesting_offline
1816
1817       line = ' '
1818
1819!
1820!--    Try to find stg package
1821       REWIND ( 11 )
1822       line = ' '
1823       DO WHILE ( INDEX( line, '&nesting_offl_parameters' ) == 0 )
1824          READ ( 11, '(A)', END=20 )  line
1825       ENDDO
1826       BACKSPACE ( 11 )
1827
1828!
1829!--    Read namelist
1830       READ ( 11, nesting_offl_parameters, ERR = 10, END = 20 )
1831
1832       GOTO 20
1833
1834 10    BACKSPACE( 11 )
1835       READ( 11 , '(A)') line
1836       CALL parin_fail_message( 'nesting_offl_parameters', line )
1837
1838 20    CONTINUE
1839
1840
1841    END SUBROUTINE nesting_offl_parin
1842
1843!------------------------------------------------------------------------------!
1844! Description:
1845! ------------
1846!> Writes information about offline nesting into HEADER file
1847!------------------------------------------------------------------------------!
1848    SUBROUTINE nesting_offl_header ( io )
1849
1850       INTEGER(iwp), INTENT(IN) ::  io !< Unit of the output file
1851
1852       WRITE ( io, 1 )
1853       IF ( nesting_offline )  THEN
1854          WRITE ( io, 3 )
1855       ELSE
1856          WRITE ( io, 2 )
1857       ENDIF
1858
18591 FORMAT (//' Offline nesting in COSMO model:'/                                &
1860              ' -------------------------------'/)
18612 FORMAT (' --> No offlince nesting is used (default) ')
18623 FORMAT (' --> Offlince nesting is used. Boundary data is read from dynamic input file ')
1863
1864    END SUBROUTINE nesting_offl_header 
1865
1866!------------------------------------------------------------------------------!
1867! Description:
1868! ------------
1869!> Allocate arrays used to read boundary data from NetCDF file and initialize
1870!> boundary data.
1871!------------------------------------------------------------------------------!
1872    SUBROUTINE nesting_offl_init
1873           
1874       INTEGER(iwp) ::  n !< running index for chemical species
1875
1876!
1877!--    Get time_utc_init from origin_date_time
1878       CALL get_date_time( 0.0_wp, second_of_day = time_utc_init )
1879
1880!--    Allocate arrays for geostrophic wind components. Arrays will
1881!--    incorporate 2 time levels in order to interpolate in between.
1882       ALLOCATE( nest_offl%ug(0:1,1:nzt) )
1883       ALLOCATE( nest_offl%vg(0:1,1:nzt) )
1884!
1885!--    Allocate arrays for reading left/right boundary values. Arrays will
1886!--    incorporate 2  time levels in order to interpolate in between. If the core has
1887!--    no boundary, allocate a dummy array, in order to enable netcdf parallel
1888!--    access. Dummy arrays will be allocated with dimension length zero.
1889       IF ( bc_dirichlet_l )  THEN
1890          ALLOCATE( nest_offl%u_left(0:1,nzb+1:nzt,nys:nyn)  )
1891          ALLOCATE( nest_offl%v_left(0:1,nzb+1:nzt,nysv:nyn) )
1892          ALLOCATE( nest_offl%w_left(0:1,nzb+1:nzt-1,nys:nyn) )
1893          IF ( humidity )       ALLOCATE( nest_offl%q_left(0:1,nzb+1:nzt,nys:nyn)  )
1894          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_left(0:1,nzb+1:nzt,nys:nyn) )
1895          IF ( air_chemistry )  ALLOCATE( nest_offl%chem_left(0:1,nzb+1:nzt,nys:nyn,&
1896                                          1:UBOUND( chem_species, 1 )) )
1897       ELSE
1898          ALLOCATE( nest_offl%u_left(1:1,1:1,1:1)  )
1899          ALLOCATE( nest_offl%v_left(1:1,1:1,1:1)  )
1900          ALLOCATE( nest_offl%w_left(1:1,1:1,1:1)  )
1901          IF ( humidity )       ALLOCATE( nest_offl%q_left(1:1,1:1,1:1)  )
1902          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_left(1:1,1:1,1:1)  )
1903          IF ( air_chemistry )  ALLOCATE( nest_offl%chem_left(1:1,1:1,1:1,     &
1904                                          1:UBOUND( chem_species, 1 )) )
1905       ENDIF
1906       IF ( bc_dirichlet_r )  THEN
1907          ALLOCATE( nest_offl%u_right(0:1,nzb+1:nzt,nys:nyn)  )
1908          ALLOCATE( nest_offl%v_right(0:1,nzb+1:nzt,nysv:nyn) )
1909          ALLOCATE( nest_offl%w_right(0:1,nzb+1:nzt-1,nys:nyn) )
1910          IF ( humidity )       ALLOCATE( nest_offl%q_right(0:1,nzb+1:nzt,nys:nyn)  )
1911          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_right(0:1,nzb+1:nzt,nys:nyn) )
1912          IF ( air_chemistry )  ALLOCATE( nest_offl%chem_right(0:1,nzb+1:nzt,nys:nyn,&
1913                                          1:UBOUND( chem_species, 1 )) )
1914       ELSE
1915          ALLOCATE( nest_offl%u_right(1:1,1:1,1:1)  )
1916          ALLOCATE( nest_offl%v_right(1:1,1:1,1:1)  )
1917          ALLOCATE( nest_offl%w_right(1:1,1:1,1:1)  )
1918          IF ( humidity )       ALLOCATE( nest_offl%q_right(1:1,1:1,1:1)  )
1919          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_right(1:1,1:1,1:1)  )
1920          IF ( air_chemistry )  ALLOCATE( nest_offl%chem_right(1:1,1:1,1:1,    &
1921                                          1:UBOUND( chem_species, 1 )) )
1922       ENDIF
1923!
1924!--    Allocate arrays for reading north/south boundary values. Arrays will
1925!--    incorporate 2  time levels in order to interpolate in between. If the core has
1926!--    no boundary, allocate a dummy array, in order to enable netcdf parallel
1927!--    access. Dummy arrays will be allocated with dimension length zero.
1928       IF ( bc_dirichlet_n )  THEN
1929          ALLOCATE( nest_offl%u_north(0:1,nzb+1:nzt,nxlu:nxr) )
1930          ALLOCATE( nest_offl%v_north(0:1,nzb+1:nzt,nxl:nxr)  )
1931          ALLOCATE( nest_offl%w_north(0:1,nzb+1:nzt-1,nxl:nxr) )
1932          IF ( humidity )       ALLOCATE( nest_offl%q_north(0:1,nzb+1:nzt,nxl:nxr)  )
1933          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_north(0:1,nzb+1:nzt,nxl:nxr) )
1934          IF ( air_chemistry )  ALLOCATE( nest_offl%chem_north(0:1,nzb+1:nzt,nxl:nxr,&
1935                                          1:UBOUND( chem_species, 1 )) )
1936       ELSE
1937          ALLOCATE( nest_offl%u_north(1:1,1:1,1:1)  )
1938          ALLOCATE( nest_offl%v_north(1:1,1:1,1:1)  )
1939          ALLOCATE( nest_offl%w_north(1:1,1:1,1:1)  )
1940          IF ( humidity )       ALLOCATE( nest_offl%q_north(1:1,1:1,1:1)  )
1941          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_north(1:1,1:1,1:1)  )
1942          IF ( air_chemistry )  ALLOCATE( nest_offl%chem_north(1:1,1:1,1:1,    &
1943                                          1:UBOUND( chem_species, 1 )) )
1944       ENDIF
1945       IF ( bc_dirichlet_s )  THEN
1946          ALLOCATE( nest_offl%u_south(0:1,nzb+1:nzt,nxlu:nxr) )
1947          ALLOCATE( nest_offl%v_south(0:1,nzb+1:nzt,nxl:nxr)  )
1948          ALLOCATE( nest_offl%w_south(0:1,nzb+1:nzt-1,nxl:nxr)    )
1949          IF ( humidity )       ALLOCATE( nest_offl%q_south(0:1,nzb+1:nzt,nxl:nxr)  )
1950          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_south(0:1,nzb+1:nzt,nxl:nxr) )
1951          IF ( air_chemistry )  ALLOCATE( nest_offl%chem_south(0:1,nzb+1:nzt,nxl:nxr,&
1952                                          1:UBOUND( chem_species, 1 )) )
1953       ELSE
1954          ALLOCATE( nest_offl%u_south(1:1,1:1,1:1)  )
1955          ALLOCATE( nest_offl%v_south(1:1,1:1,1:1)  )
1956          ALLOCATE( nest_offl%w_south(1:1,1:1,1:1)  )
1957          IF ( humidity )       ALLOCATE( nest_offl%q_south(1:1,1:1,1:1)  )
1958          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_south(1:1,1:1,1:1)  )
1959          IF ( air_chemistry )  ALLOCATE( nest_offl%chem_south(1:1,1:1,1:1,    &
1960                                          1:UBOUND( chem_species, 1 )) )
1961       ENDIF
1962!
1963!--    Allocate arrays for reading data at the top boundary. In contrast to the
1964!--    lateral boundaries, every core reads these data so that no dummy
1965!--    arrays need to be allocated.
1966       ALLOCATE( nest_offl%u_top(0:1,nys:nyn,nxlu:nxr) )
1967       ALLOCATE( nest_offl%v_top(0:1,nysv:nyn,nxl:nxr) )
1968       ALLOCATE( nest_offl%w_top(0:1,nys:nyn,nxl:nxr)  )
1969       IF ( humidity )       ALLOCATE( nest_offl%q_top(0:1,nys:nyn,nxl:nxr)  )
1970       IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_top(0:1,nys:nyn,nxl:nxr) )
1971       IF ( air_chemistry )  ALLOCATE( nest_offl%chem_top(0:1,nys:nyn,nxl:nxr, &
1972                                       1:UBOUND( chem_species, 1 )) )
1973!
1974!--    For chemical species, create the names of the variables. This is necessary
1975!--    to identify the respective variable and write it onto the correct array
1976!--    in the chem_species datatype.
1977       IF ( air_chemistry )  THEN
1978          ALLOCATE( nest_offl%chem_from_file_l(1:UBOUND( chem_species, 1 )) )
1979          ALLOCATE( nest_offl%chem_from_file_n(1:UBOUND( chem_species, 1 )) )
1980          ALLOCATE( nest_offl%chem_from_file_r(1:UBOUND( chem_species, 1 )) )
1981          ALLOCATE( nest_offl%chem_from_file_s(1:UBOUND( chem_species, 1 )) )
1982          ALLOCATE( nest_offl%chem_from_file_t(1:UBOUND( chem_species, 1 )) )
1983         
1984          ALLOCATE( nest_offl%var_names_chem_l(1:UBOUND( chem_species, 1 )) )
1985          ALLOCATE( nest_offl%var_names_chem_n(1:UBOUND( chem_species, 1 )) )
1986          ALLOCATE( nest_offl%var_names_chem_r(1:UBOUND( chem_species, 1 )) )
1987          ALLOCATE( nest_offl%var_names_chem_s(1:UBOUND( chem_species, 1 )) )
1988          ALLOCATE( nest_offl%var_names_chem_t(1:UBOUND( chem_species, 1 )) )
1989!
1990!--       Initialize flags that indicate whether the variable is on file or
1991!--       not. Please note, this is only necessary for chemistry variables.
1992          nest_offl%chem_from_file_l(:) = .FALSE.
1993          nest_offl%chem_from_file_n(:) = .FALSE.
1994          nest_offl%chem_from_file_r(:) = .FALSE.
1995          nest_offl%chem_from_file_s(:) = .FALSE.
1996          nest_offl%chem_from_file_t(:) = .FALSE.
1997         
1998          DO  n = 1, UBOUND( chem_species, 1 )
1999             nest_offl%var_names_chem_l(n) = nest_offl%char_l //               &
2000                                                  TRIM(chem_species(n)%name)
2001             nest_offl%var_names_chem_n(n) = nest_offl%char_n //               &
2002                                                  TRIM(chem_species(n)%name)
2003             nest_offl%var_names_chem_r(n) = nest_offl%char_r //               &
2004                                                  TRIM(chem_species(n)%name)
2005             nest_offl%var_names_chem_s(n) = nest_offl%char_s //               &
2006                                                  TRIM(chem_species(n)%name)
2007             nest_offl%var_names_chem_t(n) = nest_offl%char_t //               &
2008                                                  TRIM(chem_species(n)%name)
2009          ENDDO
2010       ENDIF
2011!
2012!--    Offline nesting for salsa
2013       IF ( salsa )  CALL salsa_nesting_offl_init
2014!
2015!--    Before initial data input is initiated, check if dynamic input file is
2016!--    present.
2017       IF ( .NOT. input_pids_dynamic )  THEN
2018          message_string = 'nesting_offline = .TRUE. requires dynamic '  //    &
2019                            'input file ' //                                   &
2020                            TRIM( input_file_dynamic ) // TRIM( coupling_char )
2021          CALL message( 'nesting_offl_init', 'PA0546', 1, 2, 0, 6, 0 )
2022       ENDIF
2023!
2024!--    Read COSMO data at lateral and top boundaries
2025       CALL nesting_offl_input
2026!
2027!--    Check if sufficient time steps are provided to cover the entire
2028!--    simulation. Note, dynamic input is only required for the 3D simulation,
2029!--    not for the soil/wall spinup. However, as the spinup time is added
2030!--    to the end_time, this must be considered here.
2031       IF ( end_time - spinup_time >                                           &
2032            nest_offl%time(nest_offl%nt-1) - time_utc_init )  THEN
2033          message_string = 'end_time of the simulation exceeds the ' //        &
2034                           'time dimension in the dynamic input file.'
2035          CALL message( 'nesting_offl_init', 'PA0183', 1, 2, 0, 6, 0 ) 
2036       ENDIF
2037
2038       IF ( nest_offl%time(0) /= time_utc_init )  THEN
2039          message_string = 'Offline nesting: time dimension must start at ' // &
2040                           ' time_utc_init.'
2041          CALL message( 'nesting_offl_init', 'PA0676', 1, 2, 0, 6, 0 )
2042       ENDIF
2043!
2044!--    Initialize boundary data. Please note, do not initialize boundaries in
2045!--    case of restart runs. This case the boundaries are already initialized
2046!--    and the boundary data from file would be on the wrong time level. 
2047       IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
2048          IF ( bc_dirichlet_l )  THEN
2049             u(nzb+1:nzt,nys:nyn,0)    = nest_offl%u_left(0,nzb+1:nzt,nys:nyn)
2050             v(nzb+1:nzt,nysv:nyn,-1)  = nest_offl%v_left(0,nzb+1:nzt,nysv:nyn) 
2051             w(nzb+1:nzt-1,nys:nyn,-1) = nest_offl%w_left(0,nzb+1:nzt-1,nys:nyn)
2052             IF ( .NOT. neutral )  pt(nzb+1:nzt,nys:nyn,-1) =                  &
2053                                      nest_offl%pt_left(0,nzb+1:nzt,nys:nyn)
2054             IF ( humidity      )  q(nzb+1:nzt,nys:nyn,-1)  =                  &
2055                                      nest_offl%q_left(0,nzb+1:nzt,nys:nyn)
2056             IF ( air_chemistry )  THEN
2057                DO  n = 1, UBOUND( chem_species, 1 )
2058                   IF( nest_offl%chem_from_file_l(n) )  THEN
2059                      chem_species(n)%conc(nzb+1:nzt,nys:nyn,-1) =             &
2060                                      nest_offl%chem_left(0,nzb+1:nzt,nys:nyn,n)
2061                   ENDIF
2062                ENDDO
2063             ENDIF
2064          ENDIF
2065          IF ( bc_dirichlet_r )  THEN
2066             u(nzb+1:nzt,nys:nyn,nxr+1)   = nest_offl%u_right(0,nzb+1:nzt,nys:nyn) 
2067             v(nzb+1:nzt,nysv:nyn,nxr+1)  = nest_offl%v_right(0,nzb+1:nzt,nysv:nyn)
2068             w(nzb+1:nzt-1,nys:nyn,nxr+1) = nest_offl%w_right(0,nzb+1:nzt-1,nys:nyn)
2069             IF ( .NOT. neutral )  pt(nzb+1:nzt,nys:nyn,nxr+1) =               &
2070                                      nest_offl%pt_right(0,nzb+1:nzt,nys:nyn)
2071             IF ( humidity      )  q(nzb+1:nzt,nys:nyn,nxr+1)  =               &
2072                                      nest_offl%q_right(0,nzb+1:nzt,nys:nyn)
2073             IF ( air_chemistry )  THEN
2074                DO  n = 1, UBOUND( chem_species, 1 )
2075                   IF( nest_offl%chem_from_file_r(n) )  THEN
2076                      chem_species(n)%conc(nzb+1:nzt,nys:nyn,nxr+1) =          &
2077                                      nest_offl%chem_right(0,nzb+1:nzt,nys:nyn,n)
2078                   ENDIF
2079                ENDDO
2080             ENDIF
2081          ENDIF
2082          IF ( bc_dirichlet_s )  THEN
2083             u(nzb+1:nzt,-1,nxlu:nxr)  = nest_offl%u_south(0,nzb+1:nzt,nxlu:nxr)
2084             v(nzb+1:nzt,0,nxl:nxr)    = nest_offl%v_south(0,nzb+1:nzt,nxl:nxr) 
2085             w(nzb+1:nzt-1,-1,nxl:nxr) = nest_offl%w_south(0,nzb+1:nzt-1,nxl:nxr) 
2086             IF ( .NOT. neutral )  pt(nzb+1:nzt,-1,nxl:nxr) =                  &
2087                                      nest_offl%pt_south(0,nzb+1:nzt,nxl:nxr)
2088             IF ( humidity      )  q(nzb+1:nzt,-1,nxl:nxr)  =                  &
2089                                      nest_offl%q_south(0,nzb+1:nzt,nxl:nxr) 
2090             IF ( air_chemistry )  THEN
2091                DO  n = 1, UBOUND( chem_species, 1 )
2092                   IF( nest_offl%chem_from_file_s(n) )  THEN
2093                      chem_species(n)%conc(nzb+1:nzt,-1,nxl:nxr) =             &
2094                                      nest_offl%chem_south(0,nzb+1:nzt,nxl:nxr,n)
2095                   ENDIF
2096                ENDDO
2097             ENDIF
2098          ENDIF
2099          IF ( bc_dirichlet_n )  THEN
2100             u(nzb+1:nzt,nyn+1,nxlu:nxr)  = nest_offl%u_north(0,nzb+1:nzt,nxlu:nxr)
2101             v(nzb+1:nzt,nyn+1,nxl:nxr)   = nest_offl%v_north(0,nzb+1:nzt,nxl:nxr) 
2102             w(nzb+1:nzt-1,nyn+1,nxl:nxr) = nest_offl%w_north(0,nzb+1:nzt-1,nxl:nxr)
2103             IF ( .NOT. neutral )  pt(nzb+1:nzt,nyn+1,nxl:nxr) =               &
2104                                      nest_offl%pt_north(0,nzb+1:nzt,nxl:nxr) 
2105             IF ( humidity      )  q(nzb+1:nzt,nyn+1,nxl:nxr)  =               &
2106                                      nest_offl%q_north(0,nzb+1:nzt,nxl:nxr)
2107             IF ( air_chemistry )  THEN
2108                DO  n = 1, UBOUND( chem_species, 1 )
2109                   IF( nest_offl%chem_from_file_n(n) )  THEN
2110                      chem_species(n)%conc(nzb+1:nzt,nyn+1,nxl:nxr) =          &
2111                                      nest_offl%chem_north(0,nzb+1:nzt,nxl:nxr,n)
2112                   ENDIF
2113                ENDDO
2114             ENDIF
2115          ENDIF
2116!         
2117!--       Initialize geostrophic wind components. Actually this is already done in
2118!--       init_3d_model when initializing_action = 'inifor', however, in speical
2119!--       case of user-defined initialization this will be done here again, in
2120!--       order to have a consistent initialization.
2121          ug(nzb+1:nzt) = nest_offl%ug(0,nzb+1:nzt)
2122          vg(nzb+1:nzt) = nest_offl%vg(0,nzb+1:nzt)
2123!         
2124!--       Set bottom and top boundary condition for geostrophic wind components
2125          ug(nzt+1) = ug(nzt)
2126          vg(nzt+1) = vg(nzt)
2127          ug(nzb)   = ug(nzb+1)
2128          vg(nzb)   = vg(nzb+1)
2129       ENDIF     
2130!
2131!--    After boundary data is initialized, mask topography at the
2132!--    boundaries for the velocity components.
2133       u = MERGE( u, 0.0_wp, BTEST( wall_flags_0, 1 ) )
2134       v = MERGE( v, 0.0_wp, BTEST( wall_flags_0, 2 ) )
2135       w = MERGE( w, 0.0_wp, BTEST( wall_flags_0, 3 ) )
2136!
2137!--    Initial calculation of the boundary layer depth from the prescribed
2138!--    boundary data. This is requiered for initialize the synthetic turbulence
2139!--    generator correctly.
2140       CALL nesting_offl_calc_zi
2141       
2142!
2143!--    After boundary data is initialized, ensure mass conservation. Not
2144!--    necessary in restart runs.
2145       IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
2146          CALL nesting_offl_mass_conservation
2147       ENDIF
2148
2149    END SUBROUTINE nesting_offl_init
2150   
2151!------------------------------------------------------------------------------!
2152! Description:
2153!------------------------------------------------------------------------------!
2154!> Interpolation function, used to interpolate boundary data in time.
2155!------------------------------------------------------------------------------!
2156    FUNCTION interpolate_in_time( var_t1, var_t2, fac  ) 
2157
2158       REAL(wp)            :: interpolate_in_time !< time-interpolated boundary value
2159       REAL(wp)            :: var_t1              !< boundary value at t1
2160       REAL(wp)            :: var_t2              !< boundary value at t2
2161       REAL(wp)            :: fac                 !< interpolation factor
2162
2163       interpolate_in_time = ( 1.0_wp - fac ) * var_t1 + fac * var_t2     
2164
2165    END FUNCTION interpolate_in_time
2166
2167
2168
2169 END MODULE nesting_offl_mod
Note: See TracBrowser for help on using the repository browser.