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

Last change on this file since 4372 was 4360, checked in by suehring, 5 years ago

Bugfix in output of time-averaged plant-canopy quanities; Output of plant-canopy data only where tall canopy is defined; land-surface model: fix wrong location strings; tests: update urban test case; all source code files: copyright update

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