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

Last change on this file since 4631 was 4582, checked in by suehring, 4 years ago

remove unused variable from last commit

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