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

Last change on this file since 4226 was 4226, checked in by suehring, 22 months ago

Offline nesting: data input modularized; time variable is defined relative to time_utc_init, so that input data is correctly mapped to actual model time; checks rephrased and new checks for the time dimension added; Netcdf input: routine to retrieve dimension length renamed

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