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

Last change on this file since 4346 was 4346, checked in by motisi, 16 months ago

Introduction of wall_flags_total_0, which currently sets bits based on static topography information used in wall_flags_static_0

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