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

Last change on this file since 3579 was 3579, checked in by suehring, 3 years ago

Bugfix in initialization of synthetic turbulence generator and calculation of turbulence scaling parameters in case of topography; checks for synthetic turbulence generator and offline nesting added

  • Property svn:keywords set to Id
File size: 46.6 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-2018 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: nesting_offl_mod.f90 3579 2018-11-29 15:32:39Z suehring $
27! Check implemented for offline nesting in child domain
28!
29! 3413 2018-10-24 10:28:44Z suehring
30! Keyword ID set
31!
32! 3404 2018-10-23 13:29:11Z suehring
33! Consider time-dependent geostrophic wind components in offline nesting
34!
35!
36! Initial Revision:
37! - separate offline nesting from large_scale_nudging_mod
38! - revise offline nesting, adjust for usage of synthetic turbulence generator
39! - adjust Rayleigh damping depending on the time-depending atmospheric
40!   conditions
41!
42!
43! Description:
44! ------------
45!> Offline nesting in larger-scale models. Boundary conditions for the simulation
46!> are read from NetCDF file and are prescribed onto the respective arrays.
47!> Further, a mass-flux correction is performed to maintain the mass balance.
48!--------------------------------------------------------------------------------!
49 MODULE nesting_offl_mod
50
51    USE arrays_3d,                                                             &
52        ONLY:  dzw, e, diss, pt, pt_init, q, q_init, s, u, u_init, ug, v,      &
53               v_init, vg, w, zu, zw         
54
55    USE control_parameters,                                                    &
56        ONLY:  bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, bc_dirichlet_s, &
57               dt_3d, dz, constant_diffusion, humidity, message_string,        &
58               nesting_offline, neutral, passive_scalar, rans_mode, rans_tke_e,&
59               time_since_reference_point, volume_flow
60               
61    USE cpulog,                                                                &
62        ONLY:  cpu_log, log_point
63
64    USE grid_variables
65
66    USE indices,                                                               &
67        ONLY:  nbgp, nx, nxl, nxlg, nxlu, nxr, nxrg, ny, nys,                  &
68               nysv, nysg, nyn, nyng, nzb, nz, nzt, wall_flags_0
69
70    USE kinds
71
72    USE netcdf_data_input_mod,                                                 &
73        ONLY:  nest_offl
74       
75    USE pegrid
76
77    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
78                                     !< bulk Richardson number of 0.25
79   
80    SAVE
81    PRIVATE
82!
83!-- Public subroutines
84    PUBLIC nesting_offl_bc, nesting_offl_check_parameters, nesting_offl_header,&
85           nesting_offl_init, nesting_offl_mass_conservation, nesting_offl_parin
86!
87!-- Public variables
88    PUBLIC zi_ribulk   
89
90    INTERFACE nesting_offl_bc
91       MODULE PROCEDURE nesting_offl_bc
92    END INTERFACE nesting_offl_bc
93   
94    INTERFACE nesting_offl_check_parameters
95       MODULE PROCEDURE nesting_offl_check_parameters
96    END INTERFACE nesting_offl_check_parameters
97   
98    INTERFACE nesting_offl_header
99       MODULE PROCEDURE nesting_offl_header
100    END INTERFACE nesting_offl_header
101   
102    INTERFACE nesting_offl_init
103       MODULE PROCEDURE nesting_offl_init
104    END INTERFACE nesting_offl_init
105           
106    INTERFACE nesting_offl_mass_conservation
107       MODULE PROCEDURE nesting_offl_mass_conservation
108    END INTERFACE nesting_offl_mass_conservation
109   
110    INTERFACE nesting_offl_parin
111       MODULE PROCEDURE nesting_offl_parin
112    END INTERFACE nesting_offl_parin
113
114 CONTAINS
115
116
117!------------------------------------------------------------------------------!
118! Description:
119! ------------
120!> In this subroutine a constant mass within the model domain is guaranteed.
121!> Larger-scale models may be based on a compressible equation system, which is
122!> not consistent with PALMs incompressible equation system. In order to avoid
123!> a decrease or increase of mass during the simulation, non-divergent flow
124!> through the lateral and top boundaries is compensated by the vertical wind
125!> component at the top boundary.
126!------------------------------------------------------------------------------!
127    SUBROUTINE nesting_offl_mass_conservation
128
129       IMPLICIT NONE
130
131       INTEGER(iwp) ::  i !< grid index in x-direction
132       INTEGER(iwp) ::  j !< grid index in y-direction
133       INTEGER(iwp) ::  k !< grid index in z-direction
134
135       REAL(wp) ::  d_area_t                        !< inverse of the total area of the horizontal model domain
136       REAL(wp) ::  w_correct                       !< vertical velocity increment required to compensate non-divergent flow through the boundaries
137       REAL(wp), DIMENSION(1:3) ::  volume_flow_l   !< local volume flow
138
139       
140       CALL  cpu_log( log_point(58), 'offline nesting', 'start' )
141       
142       volume_flow   = 0.0_wp
143       volume_flow_l = 0.0_wp
144
145       d_area_t = 1.0_wp / ( ( nx + 1 ) * dx * ( ny + 1 ) * dy )
146
147       IF ( bc_dirichlet_l )  THEN
148          i = nxl
149          DO  j = nys, nyn
150             DO  k = nzb+1, nzt
151                volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k) * dy   &
152                                     * MERGE( 1.0_wp, 0.0_wp,                  &
153                                              BTEST( wall_flags_0(k,j,i), 1 ) )
154             ENDDO
155          ENDDO
156       ENDIF
157       IF ( bc_dirichlet_r )  THEN
158          i = nxr+1
159          DO  j = nys, nyn
160             DO  k = nzb+1, nzt
161                volume_flow_l(1) = volume_flow_l(1) - u(k,j,i) * dzw(k) * dy   &
162                                     * MERGE( 1.0_wp, 0.0_wp,                  &
163                                              BTEST( wall_flags_0(k,j,i), 1 ) )
164             ENDDO
165          ENDDO
166       ENDIF
167       IF ( bc_dirichlet_s )  THEN
168          j = nys
169          DO  i = nxl, nxr
170             DO  k = nzb+1, nzt
171                volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k) * dx   &
172                                     * MERGE( 1.0_wp, 0.0_wp,                  &
173                                              BTEST( wall_flags_0(k,j,i), 2 ) )
174             ENDDO
175          ENDDO
176       ENDIF
177       IF ( bc_dirichlet_n )  THEN
178          j = nyn+1
179          DO  i = nxl, nxr
180             DO  k = nzb+1, nzt
181                volume_flow_l(2) = volume_flow_l(2) - v(k,j,i) * dzw(k) * dx   &
182                                     * MERGE( 1.0_wp, 0.0_wp,                  &
183                                              BTEST( wall_flags_0(k,j,i), 2 ) )
184             ENDDO
185          ENDDO
186       ENDIF
187!
188!--    Top boundary
189       k = nzt
190       DO  i = nxl, nxr
191          DO  j = nys, nyn
192             volume_flow_l(3) = volume_flow_l(3) - w(k,j,i) * dx * dy
193          ENDDO
194       ENDDO
195
196#if defined( __parallel )
197       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
198       CALL MPI_ALLREDUCE( volume_flow_l, volume_flow, 3, MPI_REAL, MPI_SUM,      &
199                           comm2d, ierr )
200#else
201       volume_flow = volume_flow_l
202#endif
203
204       w_correct = SUM( volume_flow ) * d_area_t
205
206       DO  i = nxl, nxr
207          DO  j = nys, nyn
208             DO  k = nzt, nzt + 1
209                w(k,j,i) = w(k,j,i) + w_correct
210             ENDDO
211          ENDDO
212       ENDDO
213       
214       CALL  cpu_log( log_point(58), 'offline nesting', 'stop' )
215
216    END SUBROUTINE nesting_offl_mass_conservation
217
218
219!------------------------------------------------------------------------------!
220! Description:
221! ------------
222!> Set the lateral and top boundary conditions in case the PALM domain is
223!> nested offline in a mesoscale model. Further, average boundary data and
224!> determine mean profiles, further used for correct damping in the sponge
225!> layer.
226!------------------------------------------------------------------------------!
227    SUBROUTINE nesting_offl_bc                     
228
229       IMPLICIT NONE
230
231       INTEGER(iwp) ::  i !< running index x-direction
232       INTEGER(iwp) ::  j !< running index y-direction
233       INTEGER(iwp) ::  k !< running index z-direction
234
235       REAL(wp) ::  fac_dt   !< interpolation factor
236       
237       REAL(wp), DIMENSION(nzb:nzt+1) ::  pt_ref   !< reference profile for potential temperature
238       REAL(wp), DIMENSION(nzb:nzt+1) ::  pt_ref_l !< reference profile for potential temperature on subdomain
239       REAL(wp), DIMENSION(nzb:nzt+1) ::  q_ref    !< reference profile for mixing ratio on subdomain
240       REAL(wp), DIMENSION(nzb:nzt+1) ::  q_ref_l  !< reference profile for mixing ratio on subdomain
241       REAL(wp), DIMENSION(nzb:nzt+1) ::  u_ref    !< reference profile for u-component on subdomain
242       REAL(wp), DIMENSION(nzb:nzt+1) ::  u_ref_l  !< reference profile for u-component on subdomain
243       REAL(wp), DIMENSION(nzb:nzt+1) ::  v_ref    !< reference profile for v-component on subdomain
244       REAL(wp), DIMENSION(nzb:nzt+1) ::  v_ref_l  !< reference profile for v-component on subdomain
245       
246       CALL  cpu_log( log_point(58), 'offline nesting', 'start' )     
247!
248!--    Set mean profiles, derived from boundary data, to zero
249       pt_ref   = 0.0_wp
250       q_ref    = 0.0_wp
251       u_ref    = 0.0_wp
252       v_ref    = 0.0_wp
253       
254       pt_ref_l = 0.0_wp
255       q_ref_l  = 0.0_wp
256       u_ref_l  = 0.0_wp
257       v_ref_l  = 0.0_wp
258!
259!--    Determine interpolation factor and limit it to 1. This is because
260!--    t+dt can slightly exceed time(tind_p) before boundary data is updated
261!--    again.
262       fac_dt = ( time_since_reference_point - nest_offl%time(nest_offl%tind)  &
263                + dt_3d ) /                                                    &
264           ( nest_offl%time(nest_offl%tind_p) - nest_offl%time(nest_offl%tind) )
265       fac_dt = MIN( 1.0_wp, fac_dt )
266!
267!--    Set boundary conditions of u-, v-, w-component, as well as q, and pt.
268!--    Note, boundary values at the left boundary: i=-1 (v,w,pt,q) and
269!--    i=0 (u), at the right boundary: i=nxr+1 (all), at the south boundary:
270!--    j=-1 (u,w,pt,q) and j=0 (v), at the north boundary: j=nyn+1 (all).
271!--    Please note, at the left (for u) and south (for v) boundary, values
272!--    for u and v are set also at i/j=-1, since these values are used in
273!--    boundary_conditions() to restore prognostic values.
274!--    Further, sum up data to calculate mean profiles from boundary data,
275!--    used for Rayleigh damping.
276       IF ( bc_dirichlet_l )  THEN
277
278          DO  j = nys, nyn
279             DO  k = nzb+1, nzt
280                u(k,j,0) = interpolate_in_time( nest_offl%u_left(0,k,j),       &
281                                                nest_offl%u_left(1,k,j),       &
282                                                fac_dt ) *                     &
283                             MERGE( 1.0_wp, 0.0_wp,                            &
284                                    BTEST( wall_flags_0(k,j,0), 1 ) )
285                u(k,j,-1) = u(k,j,0)
286             ENDDO
287             u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,j,0)
288          ENDDO
289
290          DO  j = nys, nyn
291             DO  k = nzb+1, nzt-1
292                w(k,j,-1) = interpolate_in_time( nest_offl%w_left(0,k,j),      &
293                                                 nest_offl%w_left(1,k,j),      &
294                                                 fac_dt ) *                    &
295                            MERGE( 1.0_wp, 0.0_wp,                             &
296                                   BTEST( wall_flags_0(k,j,-1), 3 ) )
297             ENDDO
298          ENDDO
299
300          DO  j = nysv, nyn
301             DO  k = nzb+1, nzt
302                v(k,j,-1) = interpolate_in_time( nest_offl%v_left(0,k,j),      &
303                                                 nest_offl%v_left(1,k,j),      &
304                                                 fac_dt ) *                    &
305                               MERGE( 1.0_wp, 0.0_wp,                          &
306                                      BTEST( wall_flags_0(k,j,-1), 2 ) )
307             ENDDO
308             v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,j,-1)
309          ENDDO
310
311          IF ( .NOT. neutral )  THEN
312             DO  j = nys, nyn
313                DO  k = nzb+1, nzt
314                   pt(k,j,-1) = interpolate_in_time( nest_offl%pt_left(0,k,j), &
315                                                     nest_offl%pt_left(1,k,j), &
316                                                     fac_dt )
317 
318                ENDDO
319                pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,j,-1)
320             ENDDO
321          ENDIF
322
323          IF ( humidity )  THEN
324             DO  j = nys, nyn
325                DO  k = nzb+1, nzt
326                   q(k,j,-1) = interpolate_in_time( nest_offl%q_left(0,k,j),   &
327                                                    nest_offl%q_left(1,k,j),   &
328                                                    fac_dt )
329 
330                ENDDO
331                q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,j,-1)
332             ENDDO
333          ENDIF
334
335       ENDIF
336
337       IF ( bc_dirichlet_r )  THEN
338
339          DO  j = nys, nyn
340             DO  k = nzb+1, nzt
341                u(k,j,nxr+1) = interpolate_in_time( nest_offl%u_right(0,k,j),  &
342                                                    nest_offl%u_right(1,k,j),  &
343                                                    fac_dt ) *                 &
344                             MERGE( 1.0_wp, 0.0_wp,                            &
345                                    BTEST( wall_flags_0(k,j,nxr+1), 1 ) )
346             ENDDO
347             u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,j,nxr+1)
348          ENDDO
349          DO  j = nys, nyn
350             DO  k = nzb+1, nzt-1
351                w(k,j,nxr+1) = interpolate_in_time( nest_offl%w_right(0,k,j),  &
352                                                    nest_offl%w_right(1,k,j),  &
353                                                    fac_dt ) *                 &
354                             MERGE( 1.0_wp, 0.0_wp,                            &
355                                    BTEST( wall_flags_0(k,j,nxr+1), 3 ) )
356             ENDDO
357          ENDDO
358
359          DO  j = nysv, nyn
360             DO  k = nzb+1, nzt
361                v(k,j,nxr+1) = interpolate_in_time( nest_offl%v_right(0,k,j),  &
362                                                    nest_offl%v_right(1,k,j),  &
363                                                    fac_dt ) *                 &
364                             MERGE( 1.0_wp, 0.0_wp,                            &
365                                    BTEST( wall_flags_0(k,j,nxr+1), 2 ) )
366             ENDDO
367             v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,j,nxr+1)
368          ENDDO
369
370          IF ( .NOT. neutral )  THEN
371             DO  j = nys, nyn
372                DO  k = nzb+1, nzt
373                   pt(k,j,nxr+1) = interpolate_in_time(                        &
374                                                  nest_offl%pt_right(0,k,j),   &
375                                                  nest_offl%pt_right(1,k,j),   &
376                                                  fac_dt )
377                ENDDO
378                pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,j,nxr+1)
379             ENDDO
380          ENDIF
381
382          IF ( humidity )  THEN
383             DO  j = nys, nyn
384                DO  k = nzb+1, nzt
385                   q(k,j,nxr+1) = interpolate_in_time(                         &
386                                                  nest_offl%q_right(0,k,j),    &
387                                                  nest_offl%q_right(1,k,j),    &
388                                                  fac_dt ) 
389 
390                ENDDO
391                q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,j,nxr+1)
392             ENDDO
393          ENDIF
394
395       ENDIF
396
397       IF ( bc_dirichlet_s )  THEN
398
399          DO  i = nxl, nxr
400             DO  k = nzb+1, nzt
401                v(k,0,i) = interpolate_in_time( nest_offl%v_south(0,k,i),      &
402                                                nest_offl%v_south(1,k,i),      &
403                                                fac_dt ) *                     &
404                           MERGE( 1.0_wp, 0.0_wp,                              &
405                                  BTEST( wall_flags_0(k,0,i), 2 ) )
406                v(k,-1,i) = v(k,0,i) 
407             ENDDO
408             v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,0,i)
409          ENDDO
410
411          DO  i = nxl, nxr
412             DO  k = nzb+1, nzt-1
413                w(k,-1,i) = interpolate_in_time( nest_offl%w_south(0,k,i),     &
414                                                 nest_offl%w_south(1,k,i),     &
415                                                 fac_dt ) *                    &
416                           MERGE( 1.0_wp, 0.0_wp,                              &
417                                  BTEST( wall_flags_0(k,-1,i), 3 ) )
418             ENDDO
419          ENDDO
420
421          DO  i = nxlu, nxr
422             DO  k = nzb+1, nzt
423                u(k,-1,i) = interpolate_in_time( nest_offl%u_south(0,k,i),     &
424                                                 nest_offl%u_south(1,k,i),     &
425                                                 fac_dt ) *                    &
426                           MERGE( 1.0_wp, 0.0_wp,                              &
427                                  BTEST( wall_flags_0(k,-1,i), 1 ) )
428             ENDDO
429             u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,-1,i)
430          ENDDO
431
432          IF ( .NOT. neutral )  THEN
433             DO  i = nxl, nxr
434                DO  k = nzb+1, nzt
435                   pt(k,-1,i) = interpolate_in_time(                           &
436                                                 nest_offl%pt_south(0,k,i),    &
437                                                 nest_offl%pt_south(1,k,i),    &
438                                                 fac_dt )
439 
440                ENDDO
441                pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,-1,i)
442             ENDDO
443          ENDIF
444
445          IF ( humidity )  THEN
446             DO  i = nxl, nxr
447                DO  k = nzb+1, nzt
448                   q(k,-1,i) = interpolate_in_time(                            &
449                                                 nest_offl%q_south(0,k,i),     &
450                                                 nest_offl%q_south(1,k,i),     &
451                                                 fac_dt )
452 
453                ENDDO
454                q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,-1,i)
455             ENDDO
456          ENDIF
457
458       ENDIF
459
460       IF ( bc_dirichlet_n )  THEN
461
462          DO  i = nxl, nxr
463             DO  k = nzb+1, nzt
464                v(k,nyn+1,i) = interpolate_in_time( nest_offl%v_north(0,k,i),  &
465                                                    nest_offl%v_north(1,k,i),  &
466                                                    fac_dt ) *                 &
467                               MERGE( 1.0_wp, 0.0_wp,                          &
468                                      BTEST( wall_flags_0(k,nyn+1,i), 2 ) )
469             ENDDO
470             v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,nyn+1,i)
471          ENDDO
472          DO  i = nxl, nxr
473             DO  k = nzb+1, nzt-1
474                w(k,nyn+1,i) = interpolate_in_time( nest_offl%w_north(0,k,i),  &
475                                                    nest_offl%w_north(1,k,i),  &
476                                                    fac_dt ) *                 &
477                               MERGE( 1.0_wp, 0.0_wp,                          &
478                                      BTEST( wall_flags_0(k,nyn+1,i), 3 ) )
479             ENDDO
480          ENDDO
481
482          DO  i = nxlu, nxr
483             DO  k = nzb+1, nzt
484                u(k,nyn+1,i) = interpolate_in_time( nest_offl%u_north(0,k,i),  &
485                                                    nest_offl%u_north(1,k,i),  &
486                                                    fac_dt ) *                 &
487                               MERGE( 1.0_wp, 0.0_wp,                          &
488                                      BTEST( wall_flags_0(k,nyn+1,i), 1 ) )
489
490             ENDDO
491             u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,nyn+1,i)
492          ENDDO
493
494          IF ( .NOT. neutral )  THEN
495             DO  i = nxl, nxr
496                DO  k = nzb+1, nzt
497                   pt(k,nyn+1,i) = interpolate_in_time(                        &
498                                                    nest_offl%pt_north(0,k,i), &
499                                                    nest_offl%pt_north(1,k,i), &
500                                                    fac_dt )
501 
502                ENDDO
503                pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,nyn+1,i)
504             ENDDO
505          ENDIF
506
507          IF ( humidity )  THEN
508             DO  i = nxl, nxr
509                DO  k = nzb+1, nzt
510                   q(k,nyn+1,i) = interpolate_in_time(                         &
511                                                    nest_offl%q_north(0,k,i),  &
512                                                    nest_offl%q_north(1,k,i),  &
513                                                    fac_dt )
514 
515                ENDDO
516                q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,nyn+1,i)
517             ENDDO
518          ENDIF
519
520       ENDIF
521!
522!--    Top boundary
523       DO  i = nxlu, nxr
524          DO  j = nys, nyn
525             u(nzt+1,j,i) = interpolate_in_time( nest_offl%u_top(0,j,i),       &
526                                                 nest_offl%u_top(1,j,i),       &
527                                                 fac_dt ) *                    &
528                           MERGE( 1.0_wp, 0.0_wp,                              &
529                                  BTEST( wall_flags_0(nzt+1,j,i), 1 ) )
530             u_ref_l(nzt+1) = u_ref_l(nzt+1) + u(nzt+1,j,i)
531          ENDDO
532       ENDDO
533
534       DO  i = nxl, nxr
535          DO  j = nysv, nyn
536             v(nzt+1,j,i) = interpolate_in_time( nest_offl%v_top(0,j,i),       &
537                                                 nest_offl%v_top(1,j,i),       &
538                                                 fac_dt ) *                    &
539                           MERGE( 1.0_wp, 0.0_wp,                              &
540                                  BTEST( wall_flags_0(nzt+1,j,i), 2 ) )
541             v_ref_l(nzt+1) = v_ref_l(nzt+1) + v(nzt+1,j,i)
542          ENDDO
543       ENDDO
544
545       DO  i = nxl, nxr
546          DO  j = nys, nyn
547             w(nzt,j,i) = interpolate_in_time( nest_offl%w_top(0,j,i),         &
548                                               nest_offl%w_top(1,j,i),         &
549                                               fac_dt ) *                      &
550                           MERGE( 1.0_wp, 0.0_wp,                              &
551                                  BTEST( wall_flags_0(nzt,j,i), 3 ) )
552             w(nzt+1,j,i) = w(nzt,j,i)
553          ENDDO
554       ENDDO
555
556
557       IF ( .NOT. neutral )  THEN
558          DO  i = nxl, nxr
559             DO  j = nys, nyn
560                pt(nzt+1,j,i) = interpolate_in_time( nest_offl%pt_top(0,j,i),  &
561                                                     nest_offl%pt_top(1,j,i),  &
562                                                     fac_dt )
563                pt_ref_l(nzt+1) = pt_ref_l(nzt+1) + pt(nzt+1,j,i)
564             ENDDO
565          ENDDO
566       ENDIF
567
568       IF ( humidity )  THEN
569          DO  i = nxl, nxr
570             DO  j = nys, nyn
571                q(nzt+1,j,i) = interpolate_in_time( nest_offl%q_top(0,j,i),    &
572                                                    nest_offl%q_top(1,j,i),    &
573                                                    fac_dt )
574                q_ref_l(nzt+1) = q_ref_l(nzt+1) + q(nzt+1,j,i)
575             ENDDO
576          ENDDO
577       ENDIF
578!
579!--    Moreover, set Neumann boundary condition for subgrid-scale TKE,
580!--    passive scalar, dissipation, and chemical species if required
581       IF ( rans_mode  .AND.  rans_tke_e )  THEN
582          IF (  bc_dirichlet_l )  diss(:,:,nxl-1) = diss(:,:,nxl)
583          IF (  bc_dirichlet_r )  diss(:,:,nxr+1) = diss(:,:,nxr)
584          IF (  bc_dirichlet_s )  diss(:,nys-1,:) = diss(:,nys,:)
585          IF (  bc_dirichlet_n )  diss(:,nyn+1,:) = diss(:,nyn,:)
586       ENDIF
587       IF ( .NOT. constant_diffusion )  THEN
588          IF (  bc_dirichlet_l )  e(:,:,nxl-1) = e(:,:,nxl)
589          IF (  bc_dirichlet_r )  e(:,:,nxr+1) = e(:,:,nxr)
590          IF (  bc_dirichlet_s )  e(:,nys-1,:) = e(:,nys,:)
591          IF (  bc_dirichlet_n )  e(:,nyn+1,:) = e(:,nyn,:)
592          e(nzt+1,:,:) = e(nzt,:,:)
593       ENDIF
594       IF ( passive_scalar )  THEN
595          IF (  bc_dirichlet_l )  s(:,:,nxl-1) = s(:,:,nxl)
596          IF (  bc_dirichlet_r )  s(:,:,nxr+1) = s(:,:,nxr)
597          IF (  bc_dirichlet_s )  s(:,nys-1,:) = s(:,nys,:)
598          IF (  bc_dirichlet_n )  s(:,nyn+1,:) = s(:,nyn,:)
599       ENDIF
600
601       CALL exchange_horiz( u, nbgp )
602       CALL exchange_horiz( v, nbgp )
603       CALL exchange_horiz( w, nbgp )
604       IF ( .NOT. neutral )  CALL exchange_horiz( pt, nbgp )
605       IF ( humidity      )  CALL exchange_horiz( q,  nbgp )
606!
607!--    In case of Rayleigh damping, where the profiles u_init, v_init
608!--    q_init and pt_init are still used, update these profiles from the
609!--    averaged boundary data.
610!--    But first, average these data.
611#if defined( __parallel )
612       CALL MPI_ALLREDUCE( u_ref_l, u_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM,     &
613                           comm2d, ierr )
614       CALL MPI_ALLREDUCE( v_ref_l, v_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM,     &
615                           comm2d, ierr )
616       IF ( humidity )  THEN
617          CALL MPI_ALLREDUCE( q_ref_l, q_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM,  &
618                              comm2d, ierr )
619       ENDIF
620       IF ( .NOT. neutral )  THEN
621          CALL MPI_ALLREDUCE( pt_ref_l, pt_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM,&
622                              comm2d, ierr )
623       ENDIF
624#else
625    u_ref  = u_ref_l
626    v_ref  = v_ref_l
627    IF ( humidity )       q_ref  = q_ref_l
628    IF ( .NOT. neutral )  pt_ref = pt_ref_l
629#endif
630!
631!-- Average data. Note, reference profiles up to nzt are derived from lateral
632!-- boundaries, at the model top it is derived from the top boundary. Thus,
633!-- number of input data is different from nzb:nzt compared to nzt+1.
634!-- Derived from lateral boundaries.
635    u_ref(nzb:nzt) = u_ref(nzb:nzt) / REAL( 2.0_wp * ( ny + 1 + nx     ),      &
636                                            KIND = wp ) 
637    v_ref(nzb:nzt) = v_ref(nzb:nzt) / REAL( 2.0_wp * ( ny   + nx + 1   ),      &
638                                            KIND = wp )
639    IF ( humidity )                                                            &
640       q_ref(nzb:nzt) = q_ref(nzb:nzt)   / REAL( 2.0_wp * ( ny + 1 + nx + 1 ), &
641                                                 KIND = wp )
642    IF ( .NOT. neutral )                                                       &
643       pt_ref(nzb:nzt) = pt_ref(nzb:nzt) / REAL( 2.0_wp * ( ny + 1 + nx + 1 ), &
644                                           KIND = wp )
645!
646!-- Derived from top boundary.   
647    u_ref(nzt+1) = u_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx     ), KIND = wp ) 
648    v_ref(nzt+1) = v_ref(nzt+1) / REAL( ( ny     ) * ( nx + 1 ), KIND = wp )
649    IF ( humidity )                                                            &
650       q_ref(nzt+1) = q_ref(nzt+1)   / REAL( ( ny + 1 ) * ( nx + 1 ), KIND = wp )
651    IF ( .NOT. neutral )                                                       &
652       pt_ref(nzt+1) = pt_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx + 1 ), KIND = wp )
653!
654!-- Write onto init profiles, which are used for damping
655    u_init = u_ref
656    v_init = v_ref
657    IF ( humidity      )  q_init  = q_ref
658    IF ( .NOT. neutral )  pt_init = pt_ref
659!
660!-- Set bottom boundary condition
661    IF ( humidity      )  q_init(nzb)  = q_init(nzb+1)
662    IF ( .NOT. neutral )  pt_init(nzb) = pt_init(nzb+1)
663!
664!-- Further, adjust Rayleigh damping height in case of time-changing conditions.
665!-- Therefore, calculate boundary-layer depth first.
666    CALL calc_zi
667    CALL adjust_sponge_layer
668   
669!
670!-- Update geostrophic wind components from dynamic input file.
671    DO  k = nzb+1, nzt
672       ug(k) = interpolate_in_time( nest_offl%ug(0,k), nest_offl%ug(1,k),      &
673                                    fac_dt )
674       vg(k) = interpolate_in_time( nest_offl%vg(0,k), nest_offl%vg(1,k),      &
675                                    fac_dt )
676    ENDDO
677    ug(nzt+1) = ug(nzt)
678    vg(nzt+1) = vg(nzt)
679   
680    CALL  cpu_log( log_point(58), 'offline nesting', 'stop' )
681
682    END SUBROUTINE nesting_offl_bc
683
684!------------------------------------------------------------------------------!
685! Description:
686!------------------------------------------------------------------------------!
687!> Calculates the boundary-layer depth from the boundary data, according to
688!> bulk-Richardson criterion.
689!------------------------------------------------------------------------------!
690    SUBROUTINE calc_zi
691       
692       USE basic_constants_and_equations_mod,                                  &
693           ONLY:  g
694       
695       USE kinds
696       
697       USE surface_mod,                                                        &
698           ONLY:  get_topography_top_index, get_topography_top_index_ji
699
700       IMPLICIT NONE
701
702       INTEGER(iwp) :: i            !< loop index in x-direction
703       INTEGER(iwp) :: j            !< loop index in y-direction
704       INTEGER(iwp) :: k            !< loop index in z-direction
705       INTEGER(iwp) :: k_surface    !< topography top index in z-direction
706   
707       REAL(wp) ::  ri_bulk                 !< bulk Richardson number
708       REAL(wp) ::  ri_bulk_crit = 0.25_wp  !< critical bulk Richardson number
709       REAL(wp) ::  topo_max                !< maximum topography level in model domain
710       REAL(wp) ::  topo_max_l              !< maximum topography level in subdomain
711       REAL(wp) ::  u_comp                  !< u-component
712       REAL(wp) ::  v_comp                  !< v-component
713       REAL(wp) ::  vpt_surface             !< near-surface virtual potential temperature
714       REAL(wp) ::  zi_l                    !< mean boundary-layer depth on subdomain   
715       REAL(wp) ::  zi_local                !< local boundary-layer depth 
716   
717       REAL(wp), DIMENSION(nzb:nzt+1) ::  vpt_col !< vertical profile of virtual potential temperature at (j,i)-grid point
718
719       
720!
721!--    Calculate mean boundary-layer height from boundary data.
722!--    Start with the left and right boundaries.
723       zi_l      = 0.0_wp
724       IF ( bc_dirichlet_l  .OR.  bc_dirichlet_r )  THEN
725!
726!--       Determine index along x. Please note, index indicates boundary
727!--       grid point for scalars.
728          i = MERGE( -1, nxr + 1, bc_dirichlet_l )
729   
730          DO  j = nys, nyn
731!
732!--          Determine topography top index at current (j,i) index
733             k_surface = get_topography_top_index_ji( j, i, 's' )
734!
735!--          Pre-compute surface virtual temperature. Therefore, use 2nd
736!--          prognostic level according to Heinze et al. (2017).
737             IF ( humidity )  THEN
738                vpt_surface = pt(k_surface+2,j,i) *                            &
739                            ( 1.0_wp + 0.61_wp * q(k_surface+2,j,i) )
740                vpt_col     = pt(:,j,i) * ( 1.0_wp + 0.61_wp * q(:,j,i) )
741             ELSE
742                vpt_surface = pt(k_surface+2,j,i)
743                vpt_col     = pt(:,j,i)
744             ENDIF
745!
746!--          Calculate local boundary layer height from bulk Richardson number,
747!--          i.e. the height where the bulk Richardson number exceeds its
748!--          critical value of 0.25 (according to Heinze et al., 2017).
749!--          Note, no interpolation of u- and v-component is made, as both
750!--          are mainly mean inflow profiles with very small spatial variation.
751             zi_local = 0.0_wp
752             DO  k = k_surface+1, nzt
753                u_comp = MERGE( u(k,j,i+1), u(k,j,i), bc_dirichlet_l )
754                v_comp = v(k,j,i)
755                ri_bulk = zu(k) * g / vpt_surface *                            &
756                          ( vpt_col(k) - vpt_surface ) /                       &
757                          ( u_comp * u_comp + v_comp * v_comp ) 
758                       
759                IF ( zi_local == 0.0_wp  .AND.  ri_bulk > ri_bulk_crit )       &
760                   zi_local = zu(k)           
761             ENDDO
762!
763!--          Assure that the minimum local boundary-layer depth is at least at
764!--          the second vertical grid level.
765             zi_l = zi_l + MAX( zi_local, zu(k_surface+2) )   
766             
767          ENDDO
768       
769       ENDIF
770!
771!--    Do the same at the north and south boundaries.
772       IF ( bc_dirichlet_s  .OR.  bc_dirichlet_n )  THEN
773       
774          j = MERGE( -1, nyn + 1, bc_dirichlet_s )
775       
776          DO  i = nxl, nxr
777             k_surface = get_topography_top_index_ji( j, i, 's' )
778 
779             IF ( humidity )  THEN
780                vpt_surface = pt(k_surface+2,j,i) *                            &
781                            ( 1.0_wp + 0.61_wp * q(k_surface+2,j,i) )
782                vpt_col     = pt(:,j,i) * ( 1.0_wp + 0.61_wp * q(:,j,i) )
783             ELSE
784                vpt_surface = pt(k_surface+2,j,i)
785                vpt_col  = pt(:,j,i)
786             ENDIF
787         
788             zi_local = 0.0_wp
789             DO  k = k_surface+1, nzt               
790                u_comp = u(k,j,i)
791                v_comp = MERGE( v(k,j+1,i), v(k,j,i), bc_dirichlet_s )
792                ri_bulk = zu(k) * g / vpt_surface *                            &
793                       ( vpt_col(k) - vpt_surface ) /                          &
794                       ( u_comp * u_comp + v_comp * v_comp ) 
795                       
796                IF ( zi_local == 0.0_wp  .AND.  ri_bulk > 0.25_wp )            &
797                   zi_local = zu(k)   
798             ENDDO
799             zi_l = zi_l + MAX( zi_local, zu(k_surface+2) )   
800         
801          ENDDO
802         
803       ENDIF
804   
805#if defined( __parallel )
806       CALL MPI_ALLREDUCE( zi_l, zi_ribulk, 1, MPI_REAL, MPI_SUM,              &
807                           comm2d, ierr )
808#else
809       zi_ribulk = zi_l
810#endif
811       zi_ribulk = zi_ribulk / REAL( 2 * nx + 2 * ny, KIND = wp )
812!
813!--    Finally, check if boundary layer depth is not below the any topography.
814!--    zi_ribulk will be used to adjust rayleigh damping height, i.e. the
815!--    lower level of the sponge layer, as well as to adjust the synthetic
816!--    turbulence generator accordingly. If Rayleigh damping would be applied
817!--    near buildings, etc., this would spoil the simulation results.
818       topo_max_l = zw(MAXVAL( get_topography_top_index( 's' )))
819       
820#if defined( __parallel )
821       CALL MPI_ALLREDUCE( topo_max_l, topo_max, 1, MPI_REAL, MPI_MAX,            &
822                           comm2d, ierr )
823#else
824       topo_max     = topo_max_l
825#endif
826
827       zi_ribulk = MAX( zi_ribulk, topo_max )
828
829    END SUBROUTINE calc_zi
830   
831   
832!------------------------------------------------------------------------------!
833! Description:
834!------------------------------------------------------------------------------!
835!> Adjust the height where the rayleigh damping starts, i.e. the lower level
836!> of the sponge layer.
837!------------------------------------------------------------------------------!
838    SUBROUTINE adjust_sponge_layer
839       
840       USE arrays_3d,                                                          &
841           ONLY:  rdf, rdf_sc, zu
842       
843       USE basic_constants_and_equations_mod,                                  &
844           ONLY:  pi
845           
846       USE control_parameters,                                                 &
847           ONLY:  rayleigh_damping_height, rayleigh_damping_factor
848       
849       USE kinds
850
851       IMPLICIT NONE
852
853       INTEGER(iwp) :: k   !< loop index in z-direction
854   
855       REAL(wp) ::  rdh    !< updated Rayleigh damping height
856 
857   
858       IF ( rayleigh_damping_height > 0.0_wp  .AND.                            &
859            rayleigh_damping_factor > 0.0_wp )  THEN
860!
861!--       Update Rayleigh-damping height and re-calculate height-depending
862!--       damping coefficients.
863!--       Assure that rayleigh damping starts well above the boundary layer.   
864          rdh = MIN( MAX( zi_ribulk * 1.3_wp, 10.0_wp * dz(1) ),               & 
865                     0.8_wp * zu(nzt), rayleigh_damping_height )
866!
867!--       Update Rayleigh damping factor
868          DO  k = nzb+1, nzt
869             IF ( zu(k) >= rdh )  THEN
870                rdf(k) = rayleigh_damping_factor *                             &
871                                          ( SIN( pi * 0.5_wp * ( zu(k) - rdh ) &
872                                        / ( zu(nzt) - rdh ) )                  &
873                                          )**2
874             ENDIF
875          ENDDO
876          rdf_sc = rdf
877       
878       ENDIF
879
880    END SUBROUTINE adjust_sponge_layer
881   
882!------------------------------------------------------------------------------!
883! Description:
884! ------------
885!> Performs consistency checks
886!------------------------------------------------------------------------------!
887    SUBROUTINE nesting_offl_check_parameters
888   
889       USE control_parameters,                                                    &
890        ONLY:  child_domain, message_string, nesting_offline
891
892       IMPLICIT NONE
893!
894!--    Perform checks
895       IF ( nesting_offline  .AND.  child_domain )  THEN
896          message_string = 'Offline nesting is only applicable in root model.'
897          CALL message( 'stg_check_parameters', 'PA0622', 1, 2, 0, 6, 0 )       
898       ENDIF
899
900
901    END SUBROUTINE nesting_offl_check_parameters
902   
903!------------------------------------------------------------------------------!
904! Description:
905! ------------
906!> Reads the parameter list nesting_offl_parameters
907!------------------------------------------------------------------------------!
908    SUBROUTINE nesting_offl_parin
909
910       IMPLICIT NONE
911       
912       CHARACTER (LEN=80) ::  line   !< dummy string that contains the current line of the parameter file
913
914
915       NAMELIST /nesting_offl_parameters/   nesting_offline
916
917       line = ' '
918
919!
920!--    Try to find stg package
921       REWIND ( 11 )
922       line = ' '
923       DO WHILE ( INDEX( line, '&nesting_offl_parameters' ) == 0 )
924          READ ( 11, '(A)', END=20 )  line
925       ENDDO
926       BACKSPACE ( 11 )
927
928!
929!--    Read namelist
930       READ ( 11, nesting_offl_parameters, ERR = 10, END = 20 )
931
932       GOTO 20
933
934 10    BACKSPACE( 11 )
935       READ( 11 , '(A)') line
936       CALL parin_fail_message( 'nesting_offl_parameters', line )
937
938 20    CONTINUE
939
940
941    END SUBROUTINE nesting_offl_parin
942
943!------------------------------------------------------------------------------!
944! Description:
945! ------------
946!> Writes information about offline nesting into HEADER file
947!------------------------------------------------------------------------------!
948    SUBROUTINE nesting_offl_header ( io )
949
950       IMPLICIT NONE
951
952       INTEGER(iwp), INTENT(IN) ::  io !< Unit of the output file
953
954       WRITE ( io, 1 )
955       IF ( nesting_offline )  THEN
956          WRITE ( io, 3 )
957       ELSE
958          WRITE ( io, 2 )
959       ENDIF
960
9611 FORMAT (//' Offline nesting in COSMO model:'/                                &
962              ' -------------------------------'/)
9632 FORMAT (' --> No offlince nesting is used (default) ')
9643 FORMAT (' --> Offlince nesting is used. Boundary data is read from dynamic input file ')
965
966    END SUBROUTINE nesting_offl_header
967
968!------------------------------------------------------------------------------!
969! Description:
970! ------------
971!> Allocate arrays used to read boundary data from NetCDF file and initialize
972!> boundary data.
973!------------------------------------------------------------------------------!
974    SUBROUTINE nesting_offl_init
975   
976       USE netcdf_data_input_mod,                                              &
977           ONLY:  netcdf_data_input_offline_nesting 
978
979       IMPLICIT NONE
980
981
982!--    Allocate arrays for geostrophic wind components. Arrays will
983!--    incorporate 2 time levels in order to interpolate in between.
984       ALLOCATE( nest_offl%ug(0:1,1:nzt) )
985       ALLOCATE( nest_offl%vg(0:1,1:nzt) )
986!
987!--    Allocate arrays for reading boundary values. Arrays will incorporate 2
988!--    time levels in order to interpolate in between.
989       IF ( bc_dirichlet_l )  THEN
990          ALLOCATE( nest_offl%u_left(0:1,nzb+1:nzt,nys:nyn)  )
991          ALLOCATE( nest_offl%v_left(0:1,nzb+1:nzt,nysv:nyn) )
992          ALLOCATE( nest_offl%w_left(0:1,nzb+1:nzt-1,nys:nyn) )
993          IF ( humidity )       ALLOCATE( nest_offl%q_left(0:1,nzb+1:nzt,nys:nyn)  )
994          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_left(0:1,nzb+1:nzt,nys:nyn) )
995       ENDIF
996       IF ( bc_dirichlet_r )  THEN
997          ALLOCATE( nest_offl%u_right(0:1,nzb+1:nzt,nys:nyn)  )
998          ALLOCATE( nest_offl%v_right(0:1,nzb+1:nzt,nysv:nyn) )
999          ALLOCATE( nest_offl%w_right(0:1,nzb+1:nzt-1,nys:nyn) )
1000          IF ( humidity )       ALLOCATE( nest_offl%q_right(0:1,nzb+1:nzt,nys:nyn)  )
1001          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_right(0:1,nzb+1:nzt,nys:nyn) )
1002       ENDIF
1003       IF ( bc_dirichlet_n )  THEN
1004          ALLOCATE( nest_offl%u_north(0:1,nzb+1:nzt,nxlu:nxr) )
1005          ALLOCATE( nest_offl%v_north(0:1,nzb+1:nzt,nxl:nxr)  )
1006          ALLOCATE( nest_offl%w_north(0:1,nzb+1:nzt-1,nxl:nxr) )
1007          IF ( humidity )       ALLOCATE( nest_offl%q_north(0:1,nzb+1:nzt,nxl:nxr)  )
1008          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_north(0:1,nzb+1:nzt,nxl:nxr) )
1009       ENDIF
1010       IF ( bc_dirichlet_s )  THEN
1011          ALLOCATE( nest_offl%u_south(0:1,nzb+1:nzt,nxlu:nxr) )
1012          ALLOCATE( nest_offl%v_south(0:1,nzb+1:nzt,nxl:nxr)  )
1013          ALLOCATE( nest_offl%w_south(0:1,nzb+1:nzt-1,nxl:nxr)    )
1014          IF ( humidity )       ALLOCATE( nest_offl%q_south(0:1,nzb+1:nzt,nxl:nxr)  )
1015          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_south(0:1,nzb+1:nzt,nxl:nxr) )
1016       ENDIF
1017       
1018       ALLOCATE( nest_offl%u_top(0:1,nys:nyn,nxlu:nxr) )
1019       ALLOCATE( nest_offl%v_top(0:1,nysv:nyn,nxl:nxr) )
1020       ALLOCATE( nest_offl%w_top(0:1,nys:nyn,nxl:nxr)  )
1021       IF ( humidity )       ALLOCATE( nest_offl%q_top(0:1,nys:nyn,nxl:nxr)  )
1022       IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_top(0:1,nys:nyn,nxl:nxr) )
1023
1024!
1025!--    Read COSMO data at lateral and top boundaries
1026       CALL netcdf_data_input_offline_nesting
1027!
1028!--    Initialize boundary data.
1029       IF ( bc_dirichlet_l )  THEN
1030          u(nzb+1:nzt,nys:nyn,0)    = nest_offl%u_left(0,nzb+1:nzt,nys:nyn)
1031          v(nzb+1:nzt,nysv:nyn,-1)  = nest_offl%v_left(0,nzb+1:nzt,nysv:nyn) 
1032          w(nzb+1:nzt-1,nys:nyn,-1) = nest_offl%w_left(0,nzb+1:nzt-1,nys:nyn)
1033          IF ( .NOT. neutral )  pt(nzb+1:nzt,nys:nyn,-1) =                     &
1034                                   nest_offl%pt_left(0,nzb+1:nzt,nys:nyn)
1035          IF ( humidity      )  q(nzb+1:nzt,nys:nyn,-1)  =                     &
1036                                   nest_offl%q_left(0,nzb+1:nzt,nys:nyn) 
1037       ENDIF
1038       IF ( bc_dirichlet_r )  THEN
1039          u(nzb+1:nzt,nys:nyn,nxr+1)   = nest_offl%u_right(0,nzb+1:nzt,nys:nyn) 
1040          v(nzb+1:nzt,nysv:nyn,nxr+1)  = nest_offl%v_right(0,nzb+1:nzt,nysv:nyn)
1041          w(nzb+1:nzt-1,nys:nyn,nxr+1) = nest_offl%w_right(0,nzb+1:nzt-1,nys:nyn)
1042          IF ( .NOT. neutral )  pt(nzb+1:nzt,nys:nyn,nxr+1) =                  &
1043                                   nest_offl%pt_right(0,nzb+1:nzt,nys:nyn)
1044          IF ( humidity      )  q(nzb+1:nzt,nys:nyn,nxr+1)  =                  &
1045                                   nest_offl%q_right(0,nzb+1:nzt,nys:nyn)
1046       ENDIF
1047       IF ( bc_dirichlet_s )  THEN
1048          u(nzb+1:nzt,-1,nxlu:nxr)  = nest_offl%u_south(0,nzb+1:nzt,nxlu:nxr)
1049          v(nzb+1:nzt,0,nxl:nxr)    = nest_offl%v_south(0,nzb+1:nzt,nxl:nxr) 
1050          w(nzb+1:nzt-1,-1,nxl:nxr) = nest_offl%w_south(0,nzb+1:nzt-1,nxl:nxr) 
1051          IF ( .NOT. neutral )  pt(nzb+1:nzt,-1,nxl:nxr) =                     &
1052                                   nest_offl%pt_south(0,nzb+1:nzt,nxl:nxr)
1053          IF ( humidity      )  q(nzb+1:nzt,-1,nxl:nxr)  =                     &
1054                                   nest_offl%q_south(0,nzb+1:nzt,nxl:nxr) 
1055       ENDIF
1056       IF ( bc_dirichlet_n )  THEN
1057          u(nzb+1:nzt,nyn+1,nxlu:nxr)  = nest_offl%u_north(0,nzb+1:nzt,nxlu:nxr)
1058          v(nzb+1:nzt,nyn+1,nxl:nxr)   = nest_offl%v_north(0,nzb+1:nzt,nxl:nxr) 
1059          w(nzb+1:nzt-1,nyn+1,nxl:nxr) = nest_offl%w_north(0,nzb+1:nzt-1,nxl:nxr)
1060          IF ( .NOT. neutral )  pt(nzb+1:nzt,nyn+1,nxl:nxr) =                  &
1061                                   nest_offl%pt_north(0,nzb+1:nzt,nxl:nxr) 
1062          IF ( humidity      )  q(nzb+1:nzt,nyn+1,nxl:nxr)  =                  &
1063                                   nest_offl%q_north(0,nzb+1:nzt,nxl:nxr)
1064       ENDIF
1065!
1066!--    Initialize geostrophic wind components. Actually this is already done in
1067!--    init_3d_model when initializing_action = 'inifor', however, in speical
1068!--    case of user-defined initialization this will be done here again, in
1069!--    order to have a consistent initialization.
1070       ug(nzb+1:nzt) = nest_offl%ug(0,nzb+1:nzt)
1071       vg(nzb+1:nzt) = nest_offl%vg(0,nzb+1:nzt)
1072!
1073!--    Set bottom and top boundary condition for geostrophic wind components
1074       ug(nzt+1) = ug(nzt)
1075       vg(nzt+1) = vg(nzt)
1076       ug(nzb)   = ug(nzb+1)
1077       vg(nzb)   = vg(nzb+1)
1078!
1079!--    Initial calculation of the boundary layer depth from the prescribed
1080!--    boundary data. This is requiered for initialize the synthetic turbulence
1081!--    generator correctly.
1082       CALL calc_zi
1083       
1084!
1085!--    After boundary data is initialized, mask topography at the
1086!--    boundaries for the velocity components.
1087       u = MERGE( u, 0.0_wp, BTEST( wall_flags_0, 1 ) )
1088       v = MERGE( v, 0.0_wp, BTEST( wall_flags_0, 2 ) )
1089       w = MERGE( w, 0.0_wp, BTEST( wall_flags_0, 3 ) )
1090       
1091       CALL calc_zi
1092       
1093!
1094!--    After boundary data is initialized, ensure mass conservation
1095       CALL nesting_offl_mass_conservation
1096
1097    END SUBROUTINE nesting_offl_init
1098   
1099!------------------------------------------------------------------------------!
1100! Description:
1101!------------------------------------------------------------------------------!
1102!> Interpolation function, used to interpolate boundary data in time.
1103!------------------------------------------------------------------------------!
1104    FUNCTION interpolate_in_time( var_t1, var_t2, fac  ) 
1105       
1106       USE kinds
1107
1108       IMPLICIT NONE
1109
1110       REAL(wp)            :: interpolate_in_time !< time-interpolated boundary value
1111       REAL(wp)            :: var_t1              !< boundary value at t1
1112       REAL(wp)            :: var_t2              !< boundary value at t2
1113       REAL(wp)            :: fac                 !< interpolation factor
1114
1115       interpolate_in_time = ( 1.0_wp - fac ) * var_t1 + fac * var_t2     
1116
1117    END FUNCTION interpolate_in_time
1118
1119
1120
1121 END MODULE nesting_offl_mod
Note: See TracBrowser for help on using the repository browser.