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

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