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

Last change on this file since 4180 was 4180, checked in by scharf, 5 years ago

removed comments in 'Former revisions' section that are older than 01.01.2019

  • Property svn:keywords set to Id
File size: 65.8 KB
Line 
1!> @file nesting_offl_mod.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2019 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: nesting_offl_mod.f90 4180 2019-08-21 14:37:54Z scharf $
27! Additional check added.
28!
29! 4168 2019-08-16 13:50:17Z suehring
30! Replace function get_topography_top_index by topo_top_ind
31!
32! 4125 2019-07-29 13:31:44Z suehring
33! In order to enable netcdf parallel access, allocate dummy arrays for the
34! lateral boundary data on cores that actually do not belong to these
35! boundaries.
36!
37! 4079 2019-07-09 18:04:41Z suehring
38! - Set boundary condition for w at nzt+1 at the lateral boundaries, even
39!   though these won't enter the numerical solution. However, due to the mass
40!   conservation these values might some up to very large values which will
41!   occur in the run-control file
42! - Bugfix in offline nesting of chemical species
43! - Do not set Neumann conditions for TKE and passive scalar
44!
45! 4022 2019-06-12 11:52:39Z suehring
46! Detection of boundary-layer depth in stable boundary layer on basis of
47! boundary data improved
48! Routine for boundary-layer depth calculation renamed and made public
49!
50! 3987 2019-05-22 09:52:13Z kanani
51! Introduce alternative switch for debug output during timestepping
52!
53! 3964 2019-05-09 09:48:32Z suehring
54! Ensure that veloctiy term in calculation of bulk Richardson number does not
55! become zero
56!
57! 3937 2019-04-29 15:09:07Z suehring
58! Set boundary conditon on upper-left and upper-south grid point for the u- and
59! v-component, respectively.
60!
61! 3891 2019-04-12 17:52:01Z suehring
62! Bugfix, do not overwrite lateral and top boundary data in case of restart
63! runs.
64!
65! 3885 2019-04-11 11:29:34Z kanani
66! Changes related to global restructuring of location messages and introduction
67! of additional debug messages
68!
69!
70! Do local data exchange for chemistry variables only when boundary data is 
71! coming from dynamic file
72!
73! 3737 2019-02-12 16:57:06Z suehring
74! Introduce mesoscale nesting for chemical species
75!
76! 3705 2019-01-29 19:56:39Z suehring
77! Formatting adjustments
78!
79! 3704 2019-01-29 19:51:41Z suehring
80! Check implemented for offline nesting in child domain
81!
82!
83! Description:
84! ------------
85!> Offline nesting in larger-scale models. Boundary conditions for the simulation
86!> are read from NetCDF file and are prescribed onto the respective arrays.
87!> Further, a mass-flux correction is performed to maintain the mass balance.
88!--------------------------------------------------------------------------------!
89 MODULE nesting_offl_mod
90
91    USE arrays_3d,                                                             &
92        ONLY:  dzw, e, diss, pt, pt_init, q, q_init, s, u, u_init, ug, v,      &
93               v_init, vg, w, zu, zw
94                 
95    USE chem_modules,                                                          &
96        ONLY:  chem_species
97
98    USE control_parameters,                                                    &
99        ONLY:  air_chemistry,                                                  & 
100               bc_dirichlet_l,                                                 & 
101               bc_dirichlet_n,                                                 &
102               bc_dirichlet_r,                                                 &
103               bc_dirichlet_s,                                                 &
104               dt_3d,                                                          &
105               dz,                                                             &
106               constant_diffusion,                                             &
107               child_domain,                                                   &
108               debug_output_timestep,                                          &
109               end_time,                                                       &
110               humidity,                                                       &
111               initializing_actions,                                           &
112               message_string,                                                 &
113               nesting_offline,                                                &
114               neutral,                                                        &
115               passive_scalar,                                                 &
116               rans_mode,                                                      &
117               rans_tke_e,                                                     &
118               rayleigh_damping_factor,                                        & 
119               rayleigh_damping_height,                                        &
120               spinup_time,                                                    &
121               time_since_reference_point,                                     &
122               volume_flow
123               
124    USE cpulog,                                                                &
125        ONLY:  cpu_log, log_point
126
127    USE grid_variables
128
129    USE indices,                                                               &
130        ONLY:  nbgp, nx, nxl, nxlg, nxlu, nxr, nxrg, ny, nys,                  &
131               nysv, nysg, nyn, nyng, nzb, nz, nzt,                            &
132               topo_top_ind,                                                   &
133               wall_flags_0
134
135    USE kinds
136
137    USE netcdf_data_input_mod,                                                 &
138        ONLY:  nest_offl
139       
140    USE pegrid
141
142    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
143                                     !< bulk Richardson number of 0.25
144   
145    SAVE
146    PRIVATE
147!
148!-- Public subroutines
149    PUBLIC nesting_offl_bc,                                                    &
150           nesting_offl_calc_zi,                                               &
151           nesting_offl_check_parameters,                                      &
152           nesting_offl_header,                                                &
153           nesting_offl_init,                                                  &
154           nesting_offl_mass_conservation,                                     &
155           nesting_offl_parin 
156!
157!-- Public variables
158    PUBLIC zi_ribulk   
159
160    INTERFACE nesting_offl_bc
161       MODULE PROCEDURE nesting_offl_bc
162    END INTERFACE nesting_offl_bc
163   
164    INTERFACE nesting_offl_calc_zi
165       MODULE PROCEDURE nesting_offl_calc_zi
166    END INTERFACE nesting_offl_calc_zi
167   
168    INTERFACE nesting_offl_check_parameters
169       MODULE PROCEDURE nesting_offl_check_parameters
170    END INTERFACE nesting_offl_check_parameters
171   
172    INTERFACE nesting_offl_header
173       MODULE PROCEDURE nesting_offl_header
174    END INTERFACE nesting_offl_header
175   
176    INTERFACE nesting_offl_init
177       MODULE PROCEDURE nesting_offl_init
178    END INTERFACE nesting_offl_init
179           
180    INTERFACE nesting_offl_mass_conservation
181       MODULE PROCEDURE nesting_offl_mass_conservation
182    END INTERFACE nesting_offl_mass_conservation
183   
184    INTERFACE nesting_offl_parin
185       MODULE PROCEDURE nesting_offl_parin
186    END INTERFACE nesting_offl_parin
187
188 CONTAINS
189
190
191!------------------------------------------------------------------------------!
192! Description:
193! ------------
194!> In this subroutine a constant mass within the model domain is guaranteed.
195!> Larger-scale models may be based on a compressible equation system, which is
196!> not consistent with PALMs incompressible equation system. In order to avoid
197!> a decrease or increase of mass during the simulation, non-divergent flow
198!> through the lateral and top boundaries is compensated by the vertical wind
199!> component at the top boundary.
200!------------------------------------------------------------------------------!
201    SUBROUTINE nesting_offl_mass_conservation
202
203       IMPLICIT NONE
204
205       INTEGER(iwp) ::  i !< grid index in x-direction
206       INTEGER(iwp) ::  j !< grid index in y-direction
207       INTEGER(iwp) ::  k !< grid index in z-direction
208
209       REAL(wp) ::  d_area_t                        !< inverse of the total area of the horizontal model domain
210       REAL(wp) ::  w_correct                       !< vertical velocity increment required to compensate non-divergent flow through the boundaries
211       REAL(wp), DIMENSION(1:3) ::  volume_flow_l   !< local volume flow
212
213
214       IF ( debug_output_timestep )  CALL debug_message( 'nesting_offl_mass_conservation', 'start' )
215
216       CALL  cpu_log( log_point(58), 'offline nesting', 'start' )
217       
218       volume_flow   = 0.0_wp
219       volume_flow_l = 0.0_wp
220
221       d_area_t = 1.0_wp / ( ( nx + 1 ) * dx * ( ny + 1 ) * dy )
222
223       IF ( bc_dirichlet_l )  THEN
224          i = nxl
225          DO  j = nys, nyn
226             DO  k = nzb+1, nzt
227                volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k) * dy   &
228                                     * MERGE( 1.0_wp, 0.0_wp,                  &
229                                              BTEST( wall_flags_0(k,j,i), 1 ) )
230             ENDDO
231          ENDDO
232       ENDIF
233       IF ( bc_dirichlet_r )  THEN
234          i = nxr+1
235          DO  j = nys, nyn
236             DO  k = nzb+1, nzt
237                volume_flow_l(1) = volume_flow_l(1) - u(k,j,i) * dzw(k) * dy   &
238                                     * MERGE( 1.0_wp, 0.0_wp,                  &
239                                              BTEST( wall_flags_0(k,j,i), 1 ) )
240             ENDDO
241          ENDDO
242       ENDIF
243       IF ( bc_dirichlet_s )  THEN
244          j = nys
245          DO  i = nxl, nxr
246             DO  k = nzb+1, nzt
247                volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k) * dx   &
248                                     * MERGE( 1.0_wp, 0.0_wp,                  &
249                                              BTEST( wall_flags_0(k,j,i), 2 ) )
250             ENDDO
251          ENDDO
252       ENDIF
253       IF ( bc_dirichlet_n )  THEN
254          j = nyn+1
255          DO  i = nxl, nxr
256             DO  k = nzb+1, nzt
257                volume_flow_l(2) = volume_flow_l(2) - v(k,j,i) * dzw(k) * dx   &
258                                     * MERGE( 1.0_wp, 0.0_wp,                  &
259                                              BTEST( wall_flags_0(k,j,i), 2 ) )
260             ENDDO
261          ENDDO
262       ENDIF
263!
264!--    Top boundary
265       k = nzt
266       DO  i = nxl, nxr
267          DO  j = nys, nyn
268             volume_flow_l(3) = volume_flow_l(3) - w(k,j,i) * dx * dy
269          ENDDO
270       ENDDO
271
272#if defined( __parallel )
273       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
274       CALL MPI_ALLREDUCE( volume_flow_l, volume_flow, 3, MPI_REAL, MPI_SUM,      &
275                           comm2d, ierr )
276#else
277       volume_flow = volume_flow_l
278#endif
279
280       w_correct = SUM( volume_flow ) * d_area_t
281
282       DO  i = nxl, nxr
283          DO  j = nys, nyn
284             DO  k = nzt, nzt + 1
285                w(k,j,i) = w(k,j,i) + w_correct
286             ENDDO
287          ENDDO
288       ENDDO
289       
290       CALL  cpu_log( log_point(58), 'offline nesting', 'stop' )
291
292       IF ( debug_output_timestep )  CALL debug_message( 'nesting_offl_mass_conservation', 'end' )
293
294    END SUBROUTINE nesting_offl_mass_conservation
295
296
297!------------------------------------------------------------------------------!
298! Description:
299! ------------
300!> Set the lateral and top boundary conditions in case the PALM domain is
301!> nested offline in a mesoscale model. Further, average boundary data and
302!> determine mean profiles, further used for correct damping in the sponge
303!> layer.
304!------------------------------------------------------------------------------!
305    SUBROUTINE nesting_offl_bc                     
306
307       IMPLICIT NONE
308
309       INTEGER(iwp) ::  i !< running index x-direction
310       INTEGER(iwp) ::  j !< running index y-direction
311       INTEGER(iwp) ::  k !< running index z-direction
312       INTEGER(iwp) ::  n !< running index for chemical species
313
314       REAL(wp) ::  fac_dt   !< interpolation factor
315       
316       REAL(wp), DIMENSION(nzb:nzt+1) ::  pt_ref   !< reference profile for potential temperature
317       REAL(wp), DIMENSION(nzb:nzt+1) ::  pt_ref_l !< reference profile for potential temperature on subdomain
318       REAL(wp), DIMENSION(nzb:nzt+1) ::  q_ref    !< reference profile for mixing ratio on subdomain
319       REAL(wp), DIMENSION(nzb:nzt+1) ::  q_ref_l  !< reference profile for mixing ratio on subdomain
320       REAL(wp), DIMENSION(nzb:nzt+1) ::  u_ref    !< reference profile for u-component on subdomain
321       REAL(wp), DIMENSION(nzb:nzt+1) ::  u_ref_l  !< reference profile for u-component on subdomain
322       REAL(wp), DIMENSION(nzb:nzt+1) ::  v_ref    !< reference profile for v-component on subdomain
323       REAL(wp), DIMENSION(nzb:nzt+1) ::  v_ref_l  !< reference profile for v-component on subdomain
324       
325
326       IF ( debug_output_timestep )  CALL debug_message( 'nesting_offl_bc', 'start' )
327
328       CALL  cpu_log( log_point(58), 'offline nesting', 'start' )     
329!
330!--    Set mean profiles, derived from boundary data, to zero
331       pt_ref   = 0.0_wp
332       q_ref    = 0.0_wp
333       u_ref    = 0.0_wp
334       v_ref    = 0.0_wp
335       
336       pt_ref_l = 0.0_wp
337       q_ref_l  = 0.0_wp
338       u_ref_l  = 0.0_wp
339       v_ref_l  = 0.0_wp
340!
341!--    Determine interpolation factor and limit it to 1. This is because
342!--    t+dt can slightly exceed time(tind_p) before boundary data is updated
343!--    again.
344       fac_dt = ( time_since_reference_point - nest_offl%time(nest_offl%tind)  &
345                + dt_3d ) /                                                    &
346           ( nest_offl%time(nest_offl%tind_p) - nest_offl%time(nest_offl%tind) )
347       fac_dt = MIN( 1.0_wp, fac_dt )
348!
349!--    Set boundary conditions of u-, v-, w-component, as well as q, and pt.
350!--    Note, boundary values at the left boundary: i=-1 (v,w,pt,q) and
351!--    i=0 (u), at the right boundary: i=nxr+1 (all), at the south boundary:
352!--    j=-1 (u,w,pt,q) and j=0 (v), at the north boundary: j=nyn+1 (all).
353!--    Please note, at the left (for u) and south (for v) boundary, values
354!--    for u and v are set also at i/j=-1, since these values are used in
355!--    boundary_conditions() to restore prognostic values.
356!--    Further, sum up data to calculate mean profiles from boundary data,
357!--    used for Rayleigh damping.
358       IF ( bc_dirichlet_l )  THEN
359
360          DO  j = nys, nyn
361             DO  k = nzb+1, nzt
362                u(k,j,0) = interpolate_in_time( nest_offl%u_left(0,k,j),       &
363                                                nest_offl%u_left(1,k,j),       &
364                                                fac_dt ) *                     &
365                             MERGE( 1.0_wp, 0.0_wp,                            &
366                                    BTEST( wall_flags_0(k,j,0), 1 ) )
367                u(k,j,-1) = u(k,j,0)
368             ENDDO
369             u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,j,0)
370          ENDDO
371
372          DO  j = nys, nyn
373             DO  k = nzb+1, nzt-1
374                w(k,j,-1) = interpolate_in_time( nest_offl%w_left(0,k,j),      &
375                                                 nest_offl%w_left(1,k,j),      &
376                                                 fac_dt ) *                    &
377                            MERGE( 1.0_wp, 0.0_wp,                             &
378                                   BTEST( wall_flags_0(k,j,-1), 3 ) )
379             ENDDO
380             w(nzt,j,-1) = w(nzt-1,j,-1)
381          ENDDO
382
383          DO  j = nysv, nyn
384             DO  k = nzb+1, nzt
385                v(k,j,-1) = interpolate_in_time( nest_offl%v_left(0,k,j),      &
386                                                 nest_offl%v_left(1,k,j),      &
387                                                 fac_dt ) *                    &
388                               MERGE( 1.0_wp, 0.0_wp,                          &
389                                      BTEST( wall_flags_0(k,j,-1), 2 ) )
390             ENDDO
391             v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,j,-1)
392          ENDDO
393
394          IF ( .NOT. neutral )  THEN
395             DO  j = nys, nyn
396                DO  k = nzb+1, nzt
397                   pt(k,j,-1) = interpolate_in_time( nest_offl%pt_left(0,k,j), &
398                                                     nest_offl%pt_left(1,k,j), &
399                                                     fac_dt )
400 
401                ENDDO
402                pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,j,-1)
403             ENDDO
404          ENDIF
405
406          IF ( humidity )  THEN
407             DO  j = nys, nyn
408                DO  k = nzb+1, nzt
409                   q(k,j,-1) = interpolate_in_time( nest_offl%q_left(0,k,j),   &
410                                                    nest_offl%q_left(1,k,j),   &
411                                                    fac_dt )
412 
413                ENDDO
414                q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,j,-1)
415             ENDDO
416          ENDIF
417         
418          IF ( air_chemistry )  THEN
419             DO  n = 1, UBOUND( chem_species, 1 )
420                IF ( nest_offl%chem_from_file_l(n) )  THEN                   
421                   DO  j = nys, nyn
422                      DO  k = nzb+1, nzt
423                         chem_species(n)%conc(k,j,-1) = interpolate_in_time(   &
424                                                  nest_offl%chem_left(0,k,j,n),&
425                                                  nest_offl%chem_left(1,k,j,n),&
426                                                  fac_dt                   )
427                      ENDDO
428                   ENDDO
429                ENDIF
430             ENDDO
431          ENDIF
432
433       ENDIF
434
435       IF ( bc_dirichlet_r )  THEN
436
437          DO  j = nys, nyn
438             DO  k = nzb+1, nzt
439                u(k,j,nxr+1) = interpolate_in_time( nest_offl%u_right(0,k,j),  &
440                                                    nest_offl%u_right(1,k,j),  &
441                                                    fac_dt ) *                 &
442                             MERGE( 1.0_wp, 0.0_wp,                            &
443                                    BTEST( wall_flags_0(k,j,nxr+1), 1 ) )
444             ENDDO
445             u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,j,nxr+1)
446          ENDDO
447          DO  j = nys, nyn
448             DO  k = nzb+1, nzt-1
449                w(k,j,nxr+1) = interpolate_in_time( nest_offl%w_right(0,k,j),  &
450                                                    nest_offl%w_right(1,k,j),  &
451                                                    fac_dt ) *                 &
452                             MERGE( 1.0_wp, 0.0_wp,                            &
453                                    BTEST( wall_flags_0(k,j,nxr+1), 3 ) )
454             ENDDO
455             w(nzt,j,nxr+1) = w(nzt-1,j,nxr+1)
456          ENDDO
457
458          DO  j = nysv, nyn
459             DO  k = nzb+1, nzt
460                v(k,j,nxr+1) = interpolate_in_time( nest_offl%v_right(0,k,j),  &
461                                                    nest_offl%v_right(1,k,j),  &
462                                                    fac_dt ) *                 &
463                             MERGE( 1.0_wp, 0.0_wp,                            &
464                                    BTEST( wall_flags_0(k,j,nxr+1), 2 ) )
465             ENDDO
466             v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,j,nxr+1)
467          ENDDO
468
469          IF ( .NOT. neutral )  THEN
470             DO  j = nys, nyn
471                DO  k = nzb+1, nzt
472                   pt(k,j,nxr+1) = interpolate_in_time(                        &
473                                                  nest_offl%pt_right(0,k,j),   &
474                                                  nest_offl%pt_right(1,k,j),   &
475                                                  fac_dt )
476                ENDDO
477                pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,j,nxr+1)
478             ENDDO
479          ENDIF
480
481          IF ( humidity )  THEN
482             DO  j = nys, nyn
483                DO  k = nzb+1, nzt
484                   q(k,j,nxr+1) = interpolate_in_time(                         &
485                                                  nest_offl%q_right(0,k,j),    &
486                                                  nest_offl%q_right(1,k,j),    &
487                                                  fac_dt ) 
488 
489                ENDDO
490                q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,j,nxr+1)
491             ENDDO
492          ENDIF
493         
494          IF ( air_chemistry )  THEN
495             DO  n = 1, UBOUND( chem_species, 1 )
496                IF ( nest_offl%chem_from_file_r(n) )  THEN 
497                   DO  j = nys, nyn
498                      DO  k = nzb+1, nzt
499                         chem_species(n)%conc(k,j,nxr+1) = interpolate_in_time(&
500                                                 nest_offl%chem_right(0,k,j,n),&
501                                                 nest_offl%chem_right(1,k,j,n),&
502                                                 fac_dt                       )
503                      ENDDO
504                   ENDDO
505                ENDIF
506             ENDDO
507          ENDIF
508
509       ENDIF
510
511       IF ( bc_dirichlet_s )  THEN
512
513          DO  i = nxl, nxr
514             DO  k = nzb+1, nzt
515                v(k,0,i) = interpolate_in_time( nest_offl%v_south(0,k,i),      &
516                                                nest_offl%v_south(1,k,i),      &
517                                                fac_dt ) *                     &
518                           MERGE( 1.0_wp, 0.0_wp,                              &
519                                  BTEST( wall_flags_0(k,0,i), 2 ) )
520                v(k,-1,i) = v(k,0,i) 
521             ENDDO
522             v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,0,i)
523          ENDDO
524
525          DO  i = nxl, nxr
526             DO  k = nzb+1, nzt-1
527                w(k,-1,i) = interpolate_in_time( nest_offl%w_south(0,k,i),     &
528                                                 nest_offl%w_south(1,k,i),     &
529                                                 fac_dt ) *                    &
530                           MERGE( 1.0_wp, 0.0_wp,                              &
531                                  BTEST( wall_flags_0(k,-1,i), 3 ) )
532             ENDDO
533             w(nzt,-1,i) = w(nzt-1,-1,i)
534          ENDDO
535
536          DO  i = nxlu, nxr
537             DO  k = nzb+1, nzt
538                u(k,-1,i) = interpolate_in_time( nest_offl%u_south(0,k,i),     &
539                                                 nest_offl%u_south(1,k,i),     &
540                                                 fac_dt ) *                    &
541                           MERGE( 1.0_wp, 0.0_wp,                              &
542                                  BTEST( wall_flags_0(k,-1,i), 1 ) )
543             ENDDO
544             u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,-1,i)
545          ENDDO
546
547          IF ( .NOT. neutral )  THEN
548             DO  i = nxl, nxr
549                DO  k = nzb+1, nzt
550                   pt(k,-1,i) = interpolate_in_time(                           &
551                                                 nest_offl%pt_south(0,k,i),    &
552                                                 nest_offl%pt_south(1,k,i),    &
553                                                 fac_dt )
554 
555                ENDDO
556                pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,-1,i)
557             ENDDO
558          ENDIF
559
560          IF ( humidity )  THEN
561             DO  i = nxl, nxr
562                DO  k = nzb+1, nzt
563                   q(k,-1,i) = interpolate_in_time(                            &
564                                                 nest_offl%q_south(0,k,i),     &
565                                                 nest_offl%q_south(1,k,i),     &
566                                                 fac_dt )
567 
568                ENDDO
569                q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,-1,i)
570             ENDDO
571          ENDIF
572         
573          IF ( air_chemistry )  THEN
574             DO  n = 1, UBOUND( chem_species, 1 )
575                IF ( nest_offl%chem_from_file_s(n) )  THEN 
576                   DO  i = nxl, nxr
577                      DO  k = nzb+1, nzt
578                         chem_species(n)%conc(k,-1,i) = interpolate_in_time(   &
579                                                 nest_offl%chem_south(0,k,i,n),&
580                                                 nest_offl%chem_south(1,k,i,n),&
581                                                 fac_dt                    )
582                      ENDDO
583                   ENDDO
584                ENDIF
585             ENDDO
586          ENDIF
587
588       ENDIF
589
590       IF ( bc_dirichlet_n )  THEN
591
592          DO  i = nxl, nxr
593             DO  k = nzb+1, nzt
594                v(k,nyn+1,i) = interpolate_in_time( nest_offl%v_north(0,k,i),  &
595                                                    nest_offl%v_north(1,k,i),  &
596                                                    fac_dt ) *                 &
597                               MERGE( 1.0_wp, 0.0_wp,                          &
598                                      BTEST( wall_flags_0(k,nyn+1,i), 2 ) )
599             ENDDO
600             v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,nyn+1,i)
601          ENDDO
602          DO  i = nxl, nxr
603             DO  k = nzb+1, nzt-1
604                w(k,nyn+1,i) = interpolate_in_time( nest_offl%w_north(0,k,i),  &
605                                                    nest_offl%w_north(1,k,i),  &
606                                                    fac_dt ) *                 &
607                               MERGE( 1.0_wp, 0.0_wp,                          &
608                                      BTEST( wall_flags_0(k,nyn+1,i), 3 ) )
609             ENDDO
610             w(nzt,nyn+1,i) = w(nzt-1,nyn+1,i)
611          ENDDO
612
613          DO  i = nxlu, nxr
614             DO  k = nzb+1, nzt
615                u(k,nyn+1,i) = interpolate_in_time( nest_offl%u_north(0,k,i),  &
616                                                    nest_offl%u_north(1,k,i),  &
617                                                    fac_dt ) *                 &
618                               MERGE( 1.0_wp, 0.0_wp,                          &
619                                      BTEST( wall_flags_0(k,nyn+1,i), 1 ) )
620
621             ENDDO
622             u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,nyn+1,i)
623          ENDDO
624
625          IF ( .NOT. neutral )  THEN
626             DO  i = nxl, nxr
627                DO  k = nzb+1, nzt
628                   pt(k,nyn+1,i) = interpolate_in_time(                        &
629                                                    nest_offl%pt_north(0,k,i), &
630                                                    nest_offl%pt_north(1,k,i), &
631                                                    fac_dt )
632 
633                ENDDO
634                pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,nyn+1,i)
635             ENDDO
636          ENDIF
637
638          IF ( humidity )  THEN
639             DO  i = nxl, nxr
640                DO  k = nzb+1, nzt
641                   q(k,nyn+1,i) = interpolate_in_time(                         &
642                                                    nest_offl%q_north(0,k,i),  &
643                                                    nest_offl%q_north(1,k,i),  &
644                                                    fac_dt )
645 
646                ENDDO
647                q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,nyn+1,i)
648             ENDDO
649          ENDIF
650         
651          IF ( air_chemistry )  THEN
652             DO  n = 1, UBOUND( chem_species, 1 )
653                IF ( nest_offl%chem_from_file_n(n) )  THEN 
654                   DO  i = nxl, nxr
655                      DO  k = nzb+1, nzt
656                         chem_species(n)%conc(k,nyn+1,i) = interpolate_in_time(&
657                                                 nest_offl%chem_north(0,k,i,n),&
658                                                 nest_offl%chem_north(1,k,i,n),&
659                                                 fac_dt                       )
660                      ENDDO
661                   ENDDO
662                ENDIF
663             ENDDO
664          ENDIF
665
666       ENDIF
667!
668!--    Top boundary
669       DO  i = nxlu, nxr
670          DO  j = nys, nyn
671             u(nzt+1,j,i) = interpolate_in_time( nest_offl%u_top(0,j,i),       &
672                                                 nest_offl%u_top(1,j,i),       &
673                                                 fac_dt ) *                    &
674                           MERGE( 1.0_wp, 0.0_wp,                              &
675                                  BTEST( wall_flags_0(nzt+1,j,i), 1 ) )
676             u_ref_l(nzt+1) = u_ref_l(nzt+1) + u(nzt+1,j,i)
677          ENDDO
678       ENDDO
679!
680!--    For left boundary set boundary condition for u-component also at top
681!--    grid point.
682!--    Note, this has no effect on the numeric solution, only for data output.
683       IF ( bc_dirichlet_l )  u(nzt+1,:,nxl) = u(nzt+1,:,nxlu)
684
685       DO  i = nxl, nxr
686          DO  j = nysv, nyn
687             v(nzt+1,j,i) = interpolate_in_time( nest_offl%v_top(0,j,i),       &
688                                                 nest_offl%v_top(1,j,i),       &
689                                                 fac_dt ) *                    &
690                           MERGE( 1.0_wp, 0.0_wp,                              &
691                                  BTEST( wall_flags_0(nzt+1,j,i), 2 ) )
692             v_ref_l(nzt+1) = v_ref_l(nzt+1) + v(nzt+1,j,i)
693          ENDDO
694       ENDDO
695!
696!--    For south boundary set boundary condition for v-component also at top
697!--    grid point.
698!--    Note, this has no effect on the numeric solution, only for data output.
699       IF ( bc_dirichlet_s )  v(nzt+1,nys,:) = v(nzt+1,nysv,:)
700
701       DO  i = nxl, nxr
702          DO  j = nys, nyn
703             w(nzt,j,i) = interpolate_in_time( nest_offl%w_top(0,j,i),         &
704                                               nest_offl%w_top(1,j,i),         &
705                                               fac_dt ) *                      &
706                           MERGE( 1.0_wp, 0.0_wp,                              &
707                                  BTEST( wall_flags_0(nzt,j,i), 3 ) )
708             w(nzt+1,j,i) = w(nzt,j,i)
709          ENDDO
710       ENDDO
711
712
713       IF ( .NOT. neutral )  THEN
714          DO  i = nxl, nxr
715             DO  j = nys, nyn
716                pt(nzt+1,j,i) = interpolate_in_time( nest_offl%pt_top(0,j,i),  &
717                                                     nest_offl%pt_top(1,j,i),  &
718                                                     fac_dt )
719                pt_ref_l(nzt+1) = pt_ref_l(nzt+1) + pt(nzt+1,j,i)
720             ENDDO
721          ENDDO
722       ENDIF
723
724       IF ( humidity )  THEN
725          DO  i = nxl, nxr
726             DO  j = nys, nyn
727                q(nzt+1,j,i) = interpolate_in_time( nest_offl%q_top(0,j,i),    &
728                                                    nest_offl%q_top(1,j,i),    &
729                                                    fac_dt )
730                q_ref_l(nzt+1) = q_ref_l(nzt+1) + q(nzt+1,j,i)
731             ENDDO
732          ENDDO
733       ENDIF
734       
735       IF ( air_chemistry )  THEN
736          DO  n = 1, UBOUND( chem_species, 1 )
737             IF ( nest_offl%chem_from_file_t(n) )  THEN 
738                DO  i = nxl, nxr
739                   DO  j = nys, nyn
740                      chem_species(n)%conc(nzt+1,j,i) = interpolate_in_time(   &
741                                              nest_offl%chem_top(0,j,i,n),   &
742                                              nest_offl%chem_top(1,j,i,n),   &
743                                              fac_dt                       )
744                   ENDDO
745                ENDDO
746             ENDIF
747          ENDDO
748       ENDIF
749!
750!--    Moreover, set Neumann boundary condition for subgrid-scale TKE,
751!--    passive scalar, dissipation, and chemical species if required
752       IF ( rans_mode  .AND.  rans_tke_e )  THEN
753          IF (  bc_dirichlet_l )  diss(:,:,nxl-1) = diss(:,:,nxl)
754          IF (  bc_dirichlet_r )  diss(:,:,nxr+1) = diss(:,:,nxr)
755          IF (  bc_dirichlet_s )  diss(:,nys-1,:) = diss(:,nys,:)
756          IF (  bc_dirichlet_n )  diss(:,nyn+1,:) = diss(:,nyn,:)
757       ENDIF
758!        IF ( .NOT. constant_diffusion )  THEN
759!           IF (  bc_dirichlet_l )  e(:,:,nxl-1) = e(:,:,nxl)
760!           IF (  bc_dirichlet_r )  e(:,:,nxr+1) = e(:,:,nxr)
761!           IF (  bc_dirichlet_s )  e(:,nys-1,:) = e(:,nys,:)
762!           IF (  bc_dirichlet_n )  e(:,nyn+1,:) = e(:,nyn,:)
763!           e(nzt+1,:,:) = e(nzt,:,:)
764!        ENDIF
765!        IF ( passive_scalar )  THEN
766!           IF (  bc_dirichlet_l )  s(:,:,nxl-1) = s(:,:,nxl)
767!           IF (  bc_dirichlet_r )  s(:,:,nxr+1) = s(:,:,nxr)
768!           IF (  bc_dirichlet_s )  s(:,nys-1,:) = s(:,nys,:)
769!           IF (  bc_dirichlet_n )  s(:,nyn+1,:) = s(:,nyn,:)
770!        ENDIF
771
772       CALL exchange_horiz( u, nbgp )
773       CALL exchange_horiz( v, nbgp )
774       CALL exchange_horiz( w, nbgp )
775       IF ( .NOT. neutral )  CALL exchange_horiz( pt, nbgp )
776       IF ( humidity      )  CALL exchange_horiz( q,  nbgp )
777       IF ( air_chemistry )  THEN
778          DO  n = 1, UBOUND( chem_species, 1 )
779!
780!--          Do local exchange only when necessary, i.e. when data is coming
781!--          from dynamic file.
782             IF ( nest_offl%chem_from_file_t(n) )                              &
783                CALL exchange_horiz( chem_species(n)%conc, nbgp )
784          ENDDO
785       ENDIF
786!
787!--    Set top boundary condition at all horizontal grid points, also at the
788!--    lateral boundary grid points.
789       w(nzt+1,:,:) = w(nzt,:,:)       
790!
791!--    In case of Rayleigh damping, where the profiles u_init, v_init
792!--    q_init and pt_init are still used, update these profiles from the
793!--    averaged boundary data.
794!--    But first, average these data.
795#if defined( __parallel )
796       CALL MPI_ALLREDUCE( u_ref_l, u_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM,     &
797                           comm2d, ierr )
798       CALL MPI_ALLREDUCE( v_ref_l, v_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM,     &
799                           comm2d, ierr )
800       IF ( humidity )  THEN
801          CALL MPI_ALLREDUCE( q_ref_l, q_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM,  &
802                              comm2d, ierr )
803       ENDIF
804       IF ( .NOT. neutral )  THEN
805          CALL MPI_ALLREDUCE( pt_ref_l, pt_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM,&
806                              comm2d, ierr )
807       ENDIF
808#else
809       u_ref  = u_ref_l
810       v_ref  = v_ref_l
811       IF ( humidity )       q_ref  = q_ref_l
812       IF ( .NOT. neutral )  pt_ref = pt_ref_l
813#endif
814!
815!--    Average data. Note, reference profiles up to nzt are derived from lateral
816!--    boundaries, at the model top it is derived from the top boundary. Thus,
817!--    number of input data is different from nzb:nzt compared to nzt+1.
818!--    Derived from lateral boundaries.
819       u_ref(nzb:nzt) = u_ref(nzb:nzt) / REAL( 2.0_wp * ( ny + 1 + nx     ),   &
820                                               KIND = wp ) 
821       v_ref(nzb:nzt) = v_ref(nzb:nzt) / REAL( 2.0_wp * ( ny   + nx + 1   ),   &
822                                               KIND = wp )
823       IF ( humidity )                                                         &
824          q_ref(nzb:nzt) = q_ref(nzb:nzt)   / REAL( 2.0_wp *                   &
825                                                          ( ny + 1 + nx + 1 ), &
826                                                    KIND = wp )
827       IF ( .NOT. neutral )                                                    &
828          pt_ref(nzb:nzt) = pt_ref(nzb:nzt) / REAL( 2.0_wp *                   &
829                                                          ( ny + 1 + nx + 1 ), &
830                                              KIND = wp )
831!
832!--    Derived from top boundary.   
833       u_ref(nzt+1) = u_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx     ), KIND = wp ) 
834       v_ref(nzt+1) = v_ref(nzt+1) / REAL( ( ny     ) * ( nx + 1 ), KIND = wp )
835       IF ( humidity )                                                         &
836          q_ref(nzt+1) = q_ref(nzt+1)   / REAL( ( ny + 1 ) * ( nx + 1 ),       &
837                                                KIND = wp )
838       IF ( .NOT. neutral )                                                    &
839          pt_ref(nzt+1) = pt_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx + 1 ),       &
840                                                KIND = wp )
841!
842!--    Write onto init profiles, which are used for damping
843       u_init = u_ref
844       v_init = v_ref
845       IF ( humidity      )  q_init  = q_ref
846       IF ( .NOT. neutral )  pt_init = pt_ref
847!
848!--    Set bottom boundary condition
849       IF ( humidity      )  q_init(nzb)  = q_init(nzb+1)
850       IF ( .NOT. neutral )  pt_init(nzb) = pt_init(nzb+1)
851!
852!--    Further, adjust Rayleigh damping height in case of time-changing conditions.
853!--    Therefore, calculate boundary-layer depth first.
854       CALL nesting_offl_calc_zi
855       CALL adjust_sponge_layer 
856   
857!
858!--    Update geostrophic wind components from dynamic input file.
859       DO  k = nzb+1, nzt
860          ug(k) = interpolate_in_time( nest_offl%ug(0,k), nest_offl%ug(1,k),   &
861                                       fac_dt )
862          vg(k) = interpolate_in_time( nest_offl%vg(0,k), nest_offl%vg(1,k),   &
863                                       fac_dt )
864       ENDDO
865       ug(nzt+1) = ug(nzt)
866       vg(nzt+1) = vg(nzt)
867   
868       CALL  cpu_log( log_point(58), 'offline nesting', 'stop' )
869
870       IF ( debug_output_timestep )  CALL debug_message( 'nesting_offl_bc', 'end' )
871
872
873    END SUBROUTINE nesting_offl_bc
874
875!------------------------------------------------------------------------------!
876! Description:
877!------------------------------------------------------------------------------!
878!> Calculates the boundary-layer depth from the boundary data, according to
879!> bulk-Richardson criterion.
880!------------------------------------------------------------------------------!
881    SUBROUTINE nesting_offl_calc_zi
882       
883       USE basic_constants_and_equations_mod,                                  &
884           ONLY:  g
885       
886       USE kinds
887
888       IMPLICIT NONE
889
890       INTEGER(iwp) :: i                            !< loop index in x-direction
891       INTEGER(iwp) :: j                            !< loop index in y-direction
892       INTEGER(iwp) :: k                            !< loop index in z-direction
893       INTEGER(iwp) :: k_max_loc                    !< index of maximum wind speed along z-direction
894       INTEGER(iwp) :: k_surface                    !< topography top index in z-direction
895       INTEGER(iwp) :: num_boundary_gp_non_cyclic   !< number of non-cyclic boundaries, used for averaging ABL depth
896       INTEGER(iwp) :: num_boundary_gp_non_cyclic_l !< number of non-cyclic boundaries, used for averaging ABL depth
897   
898       REAL(wp) ::  ri_bulk                 !< bulk Richardson number
899       REAL(wp) ::  ri_bulk_crit = 0.25_wp  !< critical bulk Richardson number
900       REAL(wp) ::  topo_max                !< maximum topography level in model domain
901       REAL(wp) ::  topo_max_l              !< maximum topography level in subdomain
902       REAL(wp) ::  vpt_surface             !< near-surface virtual potential temperature
903       REAL(wp) ::  zi_l                    !< mean boundary-layer depth on subdomain   
904       REAL(wp) ::  zi_local                !< local boundary-layer depth 
905   
906       REAL(wp), DIMENSION(nzb:nzt+1) ::  vpt_col !< vertical profile of virtual potential temperature at (j,i)-grid point
907       REAL(wp), DIMENSION(nzb:nzt+1) ::  uv_abs  !< vertical profile of horizontal wind speed at (j,i)-grid point
908
909       
910!
911!--    Calculate mean boundary-layer height from boundary data.
912!--    Start with the left and right boundaries.
913       zi_l      = 0.0_wp
914       num_boundary_gp_non_cyclic_l = 0
915       IF ( bc_dirichlet_l  .OR.  bc_dirichlet_r )  THEN
916!
917!--       Sum-up and store number of boundary grid points used for averaging
918!--       ABL depth
919          num_boundary_gp_non_cyclic_l = num_boundary_gp_non_cyclic_l +        &
920                                         nxr - nxl + 1
921!
922!--       Determine index along x. Please note, index indicates boundary
923!--       grid point for scalars.
924          i = MERGE( -1, nxr + 1, bc_dirichlet_l )
925   
926          DO  j = nys, nyn
927!
928!--          Determine topography top index at current (j,i) index
929             k_surface = topo_top_ind(j,i,0)
930!
931!--          Pre-compute surface virtual temperature. Therefore, use 2nd
932!--          prognostic level according to Heinze et al. (2017).
933             IF ( humidity )  THEN
934                vpt_surface = pt(k_surface+2,j,i) *                            &
935                            ( 1.0_wp + 0.61_wp * q(k_surface+2,j,i) )
936                vpt_col     = pt(:,j,i) * ( 1.0_wp + 0.61_wp * q(:,j,i) )
937             ELSE
938                vpt_surface = pt(k_surface+2,j,i)
939                vpt_col     = pt(:,j,i)
940             ENDIF
941!
942!--          Calculate local boundary layer height from bulk Richardson number,
943!--          i.e. the height where the bulk Richardson number exceeds its
944!--          critical value of 0.25 (according to Heinze et al., 2017).
945!--          Note, no interpolation of u- and v-component is made, as both
946!--          are mainly mean inflow profiles with very small spatial variation.
947!--          Add a safety factor in case the velocity term becomes zero. This
948!--          may happen if overhanging 3D structures are directly located at
949!--          the boundary, where velocity inside the building is zero
950!--          (k_surface is the index of the lowest upward-facing surface).
951             uv_abs(:) = SQRT( MERGE( u(:,j,i+1), u(:,j,i),                   &
952                                      bc_dirichlet_l )**2 +                   &
953                               v(:,j,i)**2 )
954!
955!--          Determine index of the maximum wind speed                               
956             k_max_loc = MAXLOC( uv_abs(:), DIM = 1 ) - 1
957
958             zi_local = 0.0_wp
959             DO  k = k_surface+1, nzt
960                ri_bulk = zu(k) * g / vpt_surface *                            &
961                          ( vpt_col(k) - vpt_surface ) /                       &
962                          ( uv_abs(k) + 1E-5_wp ) 
963!
964!--             Check if critical Richardson number is exceeded. Further, check
965!--             if there is a maxium in the wind profile in order to detect also
966!--             ABL heights in the stable boundary layer.
967                IF ( zi_local == 0.0_wp  .AND.                                 &
968                     ( ri_bulk > ri_bulk_crit  .OR.  k == k_max_loc ) )        &
969                   zi_local = zu(k)           
970             ENDDO
971!
972!--          Assure that the minimum local boundary-layer depth is at least at
973!--          the second vertical grid level.
974             zi_l = zi_l + MAX( zi_local, zu(k_surface+2) )   
975             
976          ENDDO
977       
978       ENDIF
979!
980!--    Do the same at the north and south boundaries.
981       IF ( bc_dirichlet_s  .OR.  bc_dirichlet_n )  THEN
982       
983          num_boundary_gp_non_cyclic_l = num_boundary_gp_non_cyclic_l +        &
984                                         nxr - nxl + 1
985       
986          j = MERGE( -1, nyn + 1, bc_dirichlet_s )
987       
988          DO  i = nxl, nxr
989             k_surface = topo_top_ind(j,i,0)
990 
991             IF ( humidity )  THEN
992                vpt_surface = pt(k_surface+2,j,i) *                            &
993                            ( 1.0_wp + 0.61_wp * q(k_surface+2,j,i) )
994                vpt_col     = pt(:,j,i) * ( 1.0_wp + 0.61_wp * q(:,j,i) )
995             ELSE
996                vpt_surface = pt(k_surface+2,j,i)
997                vpt_col  = pt(:,j,i)
998             ENDIF
999         
1000             uv_abs(:) = SQRT( u(:,j,i)**2 +                                   &
1001                               MERGE( v(:,j+1,i), v(:,j,i),                    &
1002                               bc_dirichlet_s )**2 )
1003!
1004!--          Determine index of the maximum wind speed
1005             k_max_loc = MAXLOC( uv_abs(:), DIM = 1 ) - 1
1006         
1007             zi_local = 0.0_wp
1008             DO  k = k_surface+1, nzt               
1009                ri_bulk = zu(k) * g / vpt_surface *                            &
1010                       ( vpt_col(k) - vpt_surface ) /                          &
1011                       ( uv_abs(k) + 1E-5_wp ) 
1012!
1013!--             Check if critical Richardson number is exceeded. Further, check
1014!--             if there is a maxium in the wind profile in order to detect also
1015!--             ABL heights in the stable boundary layer.                       
1016                IF ( zi_local == 0.0_wp  .AND.                                 &
1017                     ( ri_bulk > ri_bulk_crit  .OR.  k == k_max_loc ) )        &
1018                   zi_local = zu(k)   
1019             ENDDO
1020             zi_l = zi_l + MAX( zi_local, zu(k_surface+2) )   
1021         
1022          ENDDO
1023         
1024       ENDIF
1025   
1026#if defined( __parallel )
1027       CALL MPI_ALLREDUCE( zi_l, zi_ribulk, 1, MPI_REAL, MPI_SUM,              &
1028                           comm2d, ierr )
1029       CALL MPI_ALLREDUCE( num_boundary_gp_non_cyclic_l,                       &
1030                           num_boundary_gp_non_cyclic,                         &
1031                           1, MPI_INTEGER, MPI_SUM, comm2d, ierr )
1032#else
1033       zi_ribulk = zi_l
1034       num_boundary_gp_non_cyclic = num_boundary_gp_non_cyclic_l
1035#endif
1036       zi_ribulk = zi_ribulk / REAL( num_boundary_gp_non_cyclic, KIND = wp )
1037!
1038!--    Finally, check if boundary layer depth is not below the any topography.
1039!--    zi_ribulk will be used to adjust rayleigh damping height, i.e. the
1040!--    lower level of the sponge layer, as well as to adjust the synthetic
1041!--    turbulence generator accordingly. If Rayleigh damping would be applied
1042!--    near buildings, etc., this would spoil the simulation results.
1043       topo_max_l = zw(MAXVAL( topo_top_ind(nys:nyn,nxl:nxr,0) ))
1044       
1045#if defined( __parallel )
1046       CALL MPI_ALLREDUCE( topo_max_l, topo_max, 1, MPI_REAL, MPI_MAX,         &
1047                           comm2d, ierr )
1048#else
1049       topo_max     = topo_max_l
1050#endif
1051!        zi_ribulk = MAX( zi_ribulk, topo_max )
1052       
1053    END SUBROUTINE nesting_offl_calc_zi
1054   
1055   
1056!------------------------------------------------------------------------------!
1057! Description:
1058!------------------------------------------------------------------------------!
1059!> Adjust the height where the rayleigh damping starts, i.e. the lower level
1060!> of the sponge layer.
1061!------------------------------------------------------------------------------!
1062    SUBROUTINE adjust_sponge_layer
1063       
1064       USE arrays_3d,                                                          &
1065           ONLY:  rdf, rdf_sc, zu
1066       
1067       USE basic_constants_and_equations_mod,                                  &
1068           ONLY:  pi
1069       
1070       USE kinds
1071
1072       IMPLICIT NONE
1073
1074       INTEGER(iwp) :: k   !< loop index in z-direction
1075   
1076       REAL(wp) ::  rdh    !< updated Rayleigh damping height
1077 
1078   
1079       IF ( rayleigh_damping_height > 0.0_wp  .AND.                            &
1080            rayleigh_damping_factor > 0.0_wp )  THEN
1081!
1082!--       Update Rayleigh-damping height and re-calculate height-depending
1083!--       damping coefficients.
1084!--       Assure that rayleigh damping starts well above the boundary layer.   
1085          rdh = MIN( MAX( zi_ribulk * 1.3_wp, 10.0_wp * dz(1) ),               & 
1086                     0.8_wp * zu(nzt), rayleigh_damping_height )
1087!
1088!--       Update Rayleigh damping factor
1089          DO  k = nzb+1, nzt
1090             IF ( zu(k) >= rdh )  THEN
1091                rdf(k) = rayleigh_damping_factor *                             &
1092                                          ( SIN( pi * 0.5_wp * ( zu(k) - rdh ) &
1093                                        / ( zu(nzt) - rdh ) )                  &
1094                                          )**2
1095             ENDIF
1096          ENDDO
1097          rdf_sc = rdf
1098       
1099       ENDIF
1100
1101    END SUBROUTINE adjust_sponge_layer
1102   
1103!------------------------------------------------------------------------------!
1104! Description:
1105! ------------
1106!> Performs consistency checks
1107!------------------------------------------------------------------------------!
1108    SUBROUTINE nesting_offl_check_parameters 
1109
1110       IMPLICIT NONE
1111!
1112!--    Check if offline nesting is applied in nested child domain.
1113       IF ( nesting_offline  .AND.  child_domain )  THEN
1114          message_string = 'Offline nesting is only applicable in root model.'
1115          CALL message( 'offline_nesting_check_parameters', 'PA0622', 1, 2, 0, 6, 0 )       
1116       ENDIF     
1117
1118    END SUBROUTINE nesting_offl_check_parameters
1119   
1120!------------------------------------------------------------------------------!
1121! Description:
1122! ------------
1123!> Reads the parameter list nesting_offl_parameters
1124!------------------------------------------------------------------------------!
1125    SUBROUTINE nesting_offl_parin 
1126
1127       IMPLICIT NONE
1128       
1129       CHARACTER (LEN=80) ::  line   !< dummy string that contains the current line of the parameter file
1130
1131
1132       NAMELIST /nesting_offl_parameters/   nesting_offline
1133
1134       line = ' '
1135
1136!
1137!--    Try to find stg package
1138       REWIND ( 11 )
1139       line = ' '
1140       DO WHILE ( INDEX( line, '&nesting_offl_parameters' ) == 0 )
1141          READ ( 11, '(A)', END=20 )  line
1142       ENDDO
1143       BACKSPACE ( 11 )
1144
1145!
1146!--    Read namelist
1147       READ ( 11, nesting_offl_parameters, ERR = 10, END = 20 )
1148
1149       GOTO 20
1150
1151 10    BACKSPACE( 11 )
1152       READ( 11 , '(A)') line
1153       CALL parin_fail_message( 'nesting_offl_parameters', line )
1154
1155 20    CONTINUE
1156
1157
1158    END SUBROUTINE nesting_offl_parin
1159
1160!------------------------------------------------------------------------------!
1161! Description:
1162! ------------
1163!> Writes information about offline nesting into HEADER file
1164!------------------------------------------------------------------------------!
1165    SUBROUTINE nesting_offl_header ( io )
1166
1167       IMPLICIT NONE
1168
1169       INTEGER(iwp), INTENT(IN) ::  io !< Unit of the output file
1170
1171       WRITE ( io, 1 )
1172       IF ( nesting_offline )  THEN
1173          WRITE ( io, 3 )
1174       ELSE
1175          WRITE ( io, 2 )
1176       ENDIF
1177
11781 FORMAT (//' Offline nesting in COSMO model:'/                                &
1179              ' -------------------------------'/)
11802 FORMAT (' --> No offlince nesting is used (default) ')
11813 FORMAT (' --> Offlince nesting is used. Boundary data is read from dynamic input file ')
1182
1183    END SUBROUTINE nesting_offl_header 
1184
1185!------------------------------------------------------------------------------!
1186! Description:
1187! ------------
1188!> Allocate arrays used to read boundary data from NetCDF file and initialize
1189!> boundary data.
1190!------------------------------------------------------------------------------!
1191    SUBROUTINE nesting_offl_init
1192   
1193       USE netcdf_data_input_mod,                                              &
1194           ONLY:  netcdf_data_input_offline_nesting 
1195
1196       IMPLICIT NONE
1197       
1198       INTEGER(iwp) ::  n !< running index for chemical species
1199
1200
1201!--    Allocate arrays for geostrophic wind components. Arrays will
1202!--    incorporate 2 time levels in order to interpolate in between.
1203       ALLOCATE( nest_offl%ug(0:1,1:nzt) )
1204       ALLOCATE( nest_offl%vg(0:1,1:nzt) )
1205!
1206!--    Allocate arrays for reading left/right boundary values. Arrays will
1207!--    incorporate 2  time levels in order to interpolate in between. If the core has
1208!--    no boundary, allocate a dummy array, in order to enable netcdf parallel
1209!--    access. Dummy arrays will be allocated with dimension length zero.
1210       IF ( bc_dirichlet_l )  THEN
1211          ALLOCATE( nest_offl%u_left(0:1,nzb+1:nzt,nys:nyn)  )
1212          ALLOCATE( nest_offl%v_left(0:1,nzb+1:nzt,nysv:nyn) )
1213          ALLOCATE( nest_offl%w_left(0:1,nzb+1:nzt-1,nys:nyn) )
1214          IF ( humidity )       ALLOCATE( nest_offl%q_left(0:1,nzb+1:nzt,nys:nyn)  )
1215          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_left(0:1,nzb+1:nzt,nys:nyn) )
1216          IF ( air_chemistry )  ALLOCATE( nest_offl%chem_left(0:1,nzb+1:nzt,nys:nyn,&
1217                                          1:UBOUND( chem_species, 1 )) )
1218       ELSE
1219          ALLOCATE( nest_offl%u_left(1:1,1:1,1:1)  )
1220          ALLOCATE( nest_offl%v_left(1:1,1:1,1:1)  )
1221          ALLOCATE( nest_offl%w_left(1:1,1:1,1:1)  )
1222          IF ( humidity )       ALLOCATE( nest_offl%q_left(1:1,1:1,1:1)  )
1223          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_left(1:1,1:1,1:1)  )
1224          IF ( air_chemistry )  ALLOCATE( nest_offl%chem_left(1:1,1:1,1:1,     &
1225                                          1:UBOUND( chem_species, 1 )) )
1226       ENDIF
1227       IF ( bc_dirichlet_r )  THEN
1228          ALLOCATE( nest_offl%u_right(0:1,nzb+1:nzt,nys:nyn)  )
1229          ALLOCATE( nest_offl%v_right(0:1,nzb+1:nzt,nysv:nyn) )
1230          ALLOCATE( nest_offl%w_right(0:1,nzb+1:nzt-1,nys:nyn) )
1231          IF ( humidity )       ALLOCATE( nest_offl%q_right(0:1,nzb+1:nzt,nys:nyn)  )
1232          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_right(0:1,nzb+1:nzt,nys:nyn) )
1233          IF ( air_chemistry )  ALLOCATE( nest_offl%chem_right(0:1,nzb+1:nzt,nys:nyn,&
1234                                          1:UBOUND( chem_species, 1 )) )
1235       ELSE
1236          ALLOCATE( nest_offl%u_right(1:1,1:1,1:1)  )
1237          ALLOCATE( nest_offl%v_right(1:1,1:1,1:1)  )
1238          ALLOCATE( nest_offl%w_right(1:1,1:1,1:1)  )
1239          IF ( humidity )       ALLOCATE( nest_offl%q_right(1:1,1:1,1:1)  )
1240          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_right(1:1,1:1,1:1)  )
1241          IF ( air_chemistry )  ALLOCATE( nest_offl%chem_right(1:1,1:1,1:1,    &
1242                                          1:UBOUND( chem_species, 1 )) )
1243       ENDIF
1244!
1245!--    Allocate arrays for reading north/south boundary values. Arrays will
1246!--    incorporate 2  time levels in order to interpolate in between. If the core has
1247!--    no boundary, allocate a dummy array, in order to enable netcdf parallel
1248!--    access. Dummy arrays will be allocated with dimension length zero.
1249       IF ( bc_dirichlet_n )  THEN
1250          ALLOCATE( nest_offl%u_north(0:1,nzb+1:nzt,nxlu:nxr) )
1251          ALLOCATE( nest_offl%v_north(0:1,nzb+1:nzt,nxl:nxr)  )
1252          ALLOCATE( nest_offl%w_north(0:1,nzb+1:nzt-1,nxl:nxr) )
1253          IF ( humidity )       ALLOCATE( nest_offl%q_north(0:1,nzb+1:nzt,nxl:nxr)  )
1254          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_north(0:1,nzb+1:nzt,nxl:nxr) )
1255          IF ( air_chemistry )  ALLOCATE( nest_offl%chem_north(0:1,nzb+1:nzt,nxl:nxr,&
1256                                          1:UBOUND( chem_species, 1 )) )
1257       ELSE
1258          ALLOCATE( nest_offl%u_north(1:1,1:1,1:1)  )
1259          ALLOCATE( nest_offl%v_north(1:1,1:1,1:1)  )
1260          ALLOCATE( nest_offl%w_north(1:1,1:1,1:1)  )
1261          IF ( humidity )       ALLOCATE( nest_offl%q_north(1:1,1:1,1:1)  )
1262          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_north(1:1,1:1,1:1)  )
1263          IF ( air_chemistry )  ALLOCATE( nest_offl%chem_north(1:1,1:1,1:1,    &
1264                                          1:UBOUND( chem_species, 1 )) )
1265       ENDIF
1266       IF ( bc_dirichlet_s )  THEN
1267          ALLOCATE( nest_offl%u_south(0:1,nzb+1:nzt,nxlu:nxr) )
1268          ALLOCATE( nest_offl%v_south(0:1,nzb+1:nzt,nxl:nxr)  )
1269          ALLOCATE( nest_offl%w_south(0:1,nzb+1:nzt-1,nxl:nxr)    )
1270          IF ( humidity )       ALLOCATE( nest_offl%q_south(0:1,nzb+1:nzt,nxl:nxr)  )
1271          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_south(0:1,nzb+1:nzt,nxl:nxr) )
1272          IF ( air_chemistry )  ALLOCATE( nest_offl%chem_south(0:1,nzb+1:nzt,nxl:nxr,&
1273                                          1:UBOUND( chem_species, 1 )) )
1274       ELSE
1275          ALLOCATE( nest_offl%u_south(1:1,1:1,1:1)  )
1276          ALLOCATE( nest_offl%v_south(1:1,1:1,1:1)  )
1277          ALLOCATE( nest_offl%w_south(1:1,1:1,1:1)  )
1278          IF ( humidity )       ALLOCATE( nest_offl%q_south(1:1,1:1,1:1)  )
1279          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_south(1:1,1:1,1:1)  )
1280          IF ( air_chemistry )  ALLOCATE( nest_offl%chem_south(1:1,1:1,1:1,    &
1281                                          1:UBOUND( chem_species, 1 )) )
1282       ENDIF
1283!
1284!--    Allocate arrays for reading data at the top boundary. In contrast to the
1285!--    lateral boundaries, every core reads these data so that no dummy
1286!--    arrays need to be allocated.
1287       ALLOCATE( nest_offl%u_top(0:1,nys:nyn,nxlu:nxr) )
1288       ALLOCATE( nest_offl%v_top(0:1,nysv:nyn,nxl:nxr) )
1289       ALLOCATE( nest_offl%w_top(0:1,nys:nyn,nxl:nxr)  )
1290       IF ( humidity )       ALLOCATE( nest_offl%q_top(0:1,nys:nyn,nxl:nxr)  )
1291       IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_top(0:1,nys:nyn,nxl:nxr) )
1292       IF ( air_chemistry )  ALLOCATE( nest_offl%chem_top(0:1,nys:nyn,nxl:nxr, &
1293                                       1:UBOUND( chem_species, 1 )) )
1294!
1295!--    For chemical species, create the names of the variables. This is necessary
1296!--    to identify the respective variable and write it onto the correct array
1297!--    in the chem_species datatype.
1298       IF ( air_chemistry )  THEN
1299          ALLOCATE( nest_offl%chem_from_file_l(1:UBOUND( chem_species, 1 )) )
1300          ALLOCATE( nest_offl%chem_from_file_n(1:UBOUND( chem_species, 1 )) )
1301          ALLOCATE( nest_offl%chem_from_file_r(1:UBOUND( chem_species, 1 )) )
1302          ALLOCATE( nest_offl%chem_from_file_s(1:UBOUND( chem_species, 1 )) )
1303          ALLOCATE( nest_offl%chem_from_file_t(1:UBOUND( chem_species, 1 )) )
1304         
1305          ALLOCATE( nest_offl%var_names_chem_l(1:UBOUND( chem_species, 1 )) )
1306          ALLOCATE( nest_offl%var_names_chem_n(1:UBOUND( chem_species, 1 )) )
1307          ALLOCATE( nest_offl%var_names_chem_r(1:UBOUND( chem_species, 1 )) )
1308          ALLOCATE( nest_offl%var_names_chem_s(1:UBOUND( chem_species, 1 )) )
1309          ALLOCATE( nest_offl%var_names_chem_t(1:UBOUND( chem_species, 1 )) )
1310!
1311!--       Initialize flags that indicate whether the variable is on file or
1312!--       not. Please note, this is only necessary for chemistry variables.
1313          nest_offl%chem_from_file_l(:) = .FALSE.
1314          nest_offl%chem_from_file_n(:) = .FALSE.
1315          nest_offl%chem_from_file_r(:) = .FALSE.
1316          nest_offl%chem_from_file_s(:) = .FALSE.
1317          nest_offl%chem_from_file_t(:) = .FALSE.
1318         
1319          DO  n = 1, UBOUND( chem_species, 1 )
1320             nest_offl%var_names_chem_l(n) = nest_offl%char_l //               &
1321                                                  TRIM(chem_species(n)%name)
1322             nest_offl%var_names_chem_n(n) = nest_offl%char_n //               &
1323                                                  TRIM(chem_species(n)%name)
1324             nest_offl%var_names_chem_r(n) = nest_offl%char_r //               &
1325                                                  TRIM(chem_species(n)%name)
1326             nest_offl%var_names_chem_s(n) = nest_offl%char_s //               &
1327                                                  TRIM(chem_species(n)%name)
1328             nest_offl%var_names_chem_t(n) = nest_offl%char_t //               &
1329                                                  TRIM(chem_species(n)%name)
1330          ENDDO
1331       ENDIF
1332!
1333!--    Read COSMO data at lateral and top boundaries
1334       CALL netcdf_data_input_offline_nesting
1335!
1336!--    Check if sufficient time steps are provided to cover the entire
1337!--    simulation. Note, dynamic input is only required for the 3D simulation,
1338!--    not for the soil/wall spinup. However, as the spinup time is added
1339!--    to the end_time, this must be considered here.
1340       IF ( end_time - spinup_time > nest_offl%time(nest_offl%nt-1) )  THEN
1341          message_string = 'end_time > provided time in offline nesting.'
1342          CALL message( 'offline_nesting_check_parameters', 'PA0183',          &
1343                        1, 2, 0, 6, 0 ) 
1344       ENDIF
1345!
1346!--    Initialize boundary data. Please note, do not initialize boundaries in
1347!--    case of restart runs. This case the boundaries are already initialized
1348!--    and the boundary data from file would be on the wrong time level. 
1349       IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1350          IF ( bc_dirichlet_l )  THEN
1351             u(nzb+1:nzt,nys:nyn,0)    = nest_offl%u_left(0,nzb+1:nzt,nys:nyn)
1352             v(nzb+1:nzt,nysv:nyn,-1)  = nest_offl%v_left(0,nzb+1:nzt,nysv:nyn) 
1353             w(nzb+1:nzt-1,nys:nyn,-1) = nest_offl%w_left(0,nzb+1:nzt-1,nys:nyn)
1354             IF ( .NOT. neutral )  pt(nzb+1:nzt,nys:nyn,-1) =                  &
1355                                      nest_offl%pt_left(0,nzb+1:nzt,nys:nyn)
1356             IF ( humidity      )  q(nzb+1:nzt,nys:nyn,-1)  =                  &
1357                                      nest_offl%q_left(0,nzb+1:nzt,nys:nyn)
1358             IF ( air_chemistry )  THEN
1359                DO  n = 1, UBOUND( chem_species, 1 )
1360                   IF( nest_offl%chem_from_file_l(n) )  THEN
1361                      chem_species(n)%conc(nzb+1:nzt,nys:nyn,-1) =             &
1362                                      nest_offl%chem_left(0,nzb+1:nzt,nys:nyn,n)
1363                   ENDIF
1364                ENDDO
1365             ENDIF
1366          ENDIF
1367          IF ( bc_dirichlet_r )  THEN
1368             u(nzb+1:nzt,nys:nyn,nxr+1)   = nest_offl%u_right(0,nzb+1:nzt,nys:nyn) 
1369             v(nzb+1:nzt,nysv:nyn,nxr+1)  = nest_offl%v_right(0,nzb+1:nzt,nysv:nyn)
1370             w(nzb+1:nzt-1,nys:nyn,nxr+1) = nest_offl%w_right(0,nzb+1:nzt-1,nys:nyn)
1371             IF ( .NOT. neutral )  pt(nzb+1:nzt,nys:nyn,nxr+1) =               &
1372                                      nest_offl%pt_right(0,nzb+1:nzt,nys:nyn)
1373             IF ( humidity      )  q(nzb+1:nzt,nys:nyn,nxr+1)  =               &
1374                                      nest_offl%q_right(0,nzb+1:nzt,nys:nyn)
1375             IF ( air_chemistry )  THEN
1376                DO  n = 1, UBOUND( chem_species, 1 )
1377                   IF( nest_offl%chem_from_file_r(n) )  THEN
1378                      chem_species(n)%conc(nzb+1:nzt,nys:nyn,nxr+1) =          &
1379                                      nest_offl%chem_right(0,nzb+1:nzt,nys:nyn,n)
1380                   ENDIF
1381                ENDDO
1382             ENDIF
1383          ENDIF
1384          IF ( bc_dirichlet_s )  THEN
1385             u(nzb+1:nzt,-1,nxlu:nxr)  = nest_offl%u_south(0,nzb+1:nzt,nxlu:nxr)
1386             v(nzb+1:nzt,0,nxl:nxr)    = nest_offl%v_south(0,nzb+1:nzt,nxl:nxr) 
1387             w(nzb+1:nzt-1,-1,nxl:nxr) = nest_offl%w_south(0,nzb+1:nzt-1,nxl:nxr) 
1388             IF ( .NOT. neutral )  pt(nzb+1:nzt,-1,nxl:nxr) =                  &
1389                                      nest_offl%pt_south(0,nzb+1:nzt,nxl:nxr)
1390             IF ( humidity      )  q(nzb+1:nzt,-1,nxl:nxr)  =                  &
1391                                      nest_offl%q_south(0,nzb+1:nzt,nxl:nxr) 
1392             IF ( air_chemistry )  THEN
1393                DO  n = 1, UBOUND( chem_species, 1 )
1394                   IF( nest_offl%chem_from_file_s(n) )  THEN
1395                      chem_species(n)%conc(nzb+1:nzt,-1,nxl:nxr) =             &
1396                                      nest_offl%chem_south(0,nzb+1:nzt,nxl:nxr,n)
1397                   ENDIF
1398                ENDDO
1399             ENDIF
1400          ENDIF
1401          IF ( bc_dirichlet_n )  THEN
1402             u(nzb+1:nzt,nyn+1,nxlu:nxr)  = nest_offl%u_north(0,nzb+1:nzt,nxlu:nxr)
1403             v(nzb+1:nzt,nyn+1,nxl:nxr)   = nest_offl%v_north(0,nzb+1:nzt,nxl:nxr) 
1404             w(nzb+1:nzt-1,nyn+1,nxl:nxr) = nest_offl%w_north(0,nzb+1:nzt-1,nxl:nxr)
1405             IF ( .NOT. neutral )  pt(nzb+1:nzt,nyn+1,nxl:nxr) =               &
1406                                      nest_offl%pt_north(0,nzb+1:nzt,nxl:nxr) 
1407             IF ( humidity      )  q(nzb+1:nzt,nyn+1,nxl:nxr)  =               &
1408                                      nest_offl%q_north(0,nzb+1:nzt,nxl:nxr)
1409             IF ( air_chemistry )  THEN
1410                DO  n = 1, UBOUND( chem_species, 1 )
1411                   IF( nest_offl%chem_from_file_n(n) )  THEN
1412                      chem_species(n)%conc(nzb+1:nzt,nyn+1,nxl:nxr) =          &
1413                                      nest_offl%chem_north(0,nzb+1:nzt,nxl:nxr,n)
1414                   ENDIF
1415                ENDDO
1416             ENDIF
1417          ENDIF
1418!         
1419!--       Initialize geostrophic wind components. Actually this is already done in
1420!--       init_3d_model when initializing_action = 'inifor', however, in speical
1421!--       case of user-defined initialization this will be done here again, in
1422!--       order to have a consistent initialization.
1423          ug(nzb+1:nzt) = nest_offl%ug(0,nzb+1:nzt)
1424          vg(nzb+1:nzt) = nest_offl%vg(0,nzb+1:nzt)
1425!         
1426!--       Set bottom and top boundary condition for geostrophic wind components
1427          ug(nzt+1) = ug(nzt)
1428          vg(nzt+1) = vg(nzt)
1429          ug(nzb)   = ug(nzb+1)
1430          vg(nzb)   = vg(nzb+1)
1431       ENDIF     
1432!
1433!--    After boundary data is initialized, mask topography at the
1434!--    boundaries for the velocity components.
1435       u = MERGE( u, 0.0_wp, BTEST( wall_flags_0, 1 ) )
1436       v = MERGE( v, 0.0_wp, BTEST( wall_flags_0, 2 ) )
1437       w = MERGE( w, 0.0_wp, BTEST( wall_flags_0, 3 ) )
1438!
1439!--    Initial calculation of the boundary layer depth from the prescribed
1440!--    boundary data. This is requiered for initialize the synthetic turbulence
1441!--    generator correctly.
1442       CALL nesting_offl_calc_zi
1443       
1444!
1445!--    After boundary data is initialized, ensure mass conservation. Not
1446!--    necessary in restart runs.
1447       IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1448          CALL nesting_offl_mass_conservation
1449       ENDIF
1450
1451    END SUBROUTINE nesting_offl_init
1452   
1453!------------------------------------------------------------------------------!
1454! Description:
1455!------------------------------------------------------------------------------!
1456!> Interpolation function, used to interpolate boundary data in time.
1457!------------------------------------------------------------------------------!
1458    FUNCTION interpolate_in_time( var_t1, var_t2, fac  ) 
1459       
1460       USE kinds
1461
1462       IMPLICIT NONE
1463
1464       REAL(wp)            :: interpolate_in_time !< time-interpolated boundary value
1465       REAL(wp)            :: var_t1              !< boundary value at t1
1466       REAL(wp)            :: var_t2              !< boundary value at t2
1467       REAL(wp)            :: fac                 !< interpolation factor
1468
1469       interpolate_in_time = ( 1.0_wp - fac ) * var_t1 + fac * var_t2     
1470
1471    END FUNCTION interpolate_in_time
1472
1473
1474
1475 END MODULE nesting_offl_mod
Note: See TracBrowser for help on using the repository browser.