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

Last change on this file since 4564 was 4564, checked in by raasch, 4 years ago

Vertical nesting method of Huq et al. (2019) removed

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