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

Last change on this file since 4286 was 4286, checked in by resler, 2 years ago

Merge branch resler into trunk

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