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

Last change on this file since 4178 was 4169, checked in by suehring, 2 years ago

Check added

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