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

Last change on this file since 4265 was 4231, checked in by suehring, 5 years ago

Bugfix in array deallocation

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