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

Last change on this file since 4581 was 4581, checked in by suehring, 11 months ago

mesoscale nesting: omit explicit pressure forcing via geostrophic wind components; chemistry: enable profile output of vertical fluxes; urban-surface: bugfix in initialization in case of cyclic_fill

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