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

Last change on this file since 4230 was 4230, checked in by suehring, 21 months ago

Several bugfixes: profile initialization of chemical species in restart runs; Runge-Kutta tendency array not initialized in chemistry model in restart runs; fixed determination of time indices for chemical emissions (introduced with commit -4227); Update chemistry profiles in offline nesting; initialize canopy resistances for greened building walls (even if green fraction is zero)

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