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

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

Revision of virtual-measurement module and data output enabled. Further, post-processing tool added to merge distributed virtually sampled data and to output it into NetCDF files.

  • Property svn:keywords set to Id
File size: 47.0 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! Formatting adjustments
23!
24! Former revisions:
25! -----------------
26! $Id: nesting_offl_mod.f90 3704 2019-01-29 19:51:41Z 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 *                   &
641                                                          ( ny + 1 + nx + 1 ), &
642                                                    KIND = wp )
643       IF ( .NOT. neutral )                                                    &
644          pt_ref(nzb:nzt) = pt_ref(nzb:nzt) / REAL( 2.0_wp *                   &
645                                                          ( ny + 1 + nx + 1 ), &
646                                              KIND = wp )
647!
648!--    Derived from top boundary.   
649       u_ref(nzt+1) = u_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx     ), KIND = wp ) 
650       v_ref(nzt+1) = v_ref(nzt+1) / REAL( ( ny     ) * ( nx + 1 ), KIND = wp )
651       IF ( humidity )                                                         &
652          q_ref(nzt+1) = q_ref(nzt+1)   / REAL( ( ny + 1 ) * ( nx + 1 ),       &
653                                                KIND = wp )
654       IF ( .NOT. neutral )                                                    &
655          pt_ref(nzt+1) = pt_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx + 1 ),       &
656                                                KIND = wp )
657!
658!--    Write onto init profiles, which are used for damping
659       u_init = u_ref
660       v_init = v_ref
661       IF ( humidity      )  q_init  = q_ref
662       IF ( .NOT. neutral )  pt_init = pt_ref
663!
664!--    Set bottom boundary condition
665       IF ( humidity      )  q_init(nzb)  = q_init(nzb+1)
666       IF ( .NOT. neutral )  pt_init(nzb) = pt_init(nzb+1)
667!
668!--    Further, adjust Rayleigh damping height in case of time-changing conditions.
669!--    Therefore, calculate boundary-layer depth first.
670       CALL calc_zi
671       CALL adjust_sponge_layer
672   
673!
674!--    Update geostrophic wind components from dynamic input file.
675       DO  k = nzb+1, nzt
676          ug(k) = interpolate_in_time( nest_offl%ug(0,k), nest_offl%ug(1,k),   &
677                                       fac_dt )
678          vg(k) = interpolate_in_time( nest_offl%vg(0,k), nest_offl%vg(1,k),   &
679                                       fac_dt )
680       ENDDO
681       ug(nzt+1) = ug(nzt)
682       vg(nzt+1) = vg(nzt)
683   
684       CALL  cpu_log( log_point(58), 'offline nesting', 'stop' )
685
686    END SUBROUTINE nesting_offl_bc
687
688!------------------------------------------------------------------------------!
689! Description:
690!------------------------------------------------------------------------------!
691!> Calculates the boundary-layer depth from the boundary data, according to
692!> bulk-Richardson criterion.
693!------------------------------------------------------------------------------!
694    SUBROUTINE calc_zi
695       
696       USE basic_constants_and_equations_mod,                                  &
697           ONLY:  g
698       
699       USE kinds
700       
701       USE surface_mod,                                                        &
702           ONLY:  get_topography_top_index, get_topography_top_index_ji
703
704       IMPLICIT NONE
705
706       INTEGER(iwp) :: i            !< loop index in x-direction
707       INTEGER(iwp) :: j            !< loop index in y-direction
708       INTEGER(iwp) :: k            !< loop index in z-direction
709       INTEGER(iwp) :: k_surface    !< topography top index in z-direction
710   
711       REAL(wp) ::  ri_bulk                 !< bulk Richardson number
712       REAL(wp) ::  ri_bulk_crit = 0.25_wp  !< critical bulk Richardson number
713       REAL(wp) ::  topo_max                !< maximum topography level in model domain
714       REAL(wp) ::  topo_max_l              !< maximum topography level in subdomain
715       REAL(wp) ::  u_comp                  !< u-component
716       REAL(wp) ::  v_comp                  !< v-component
717       REAL(wp) ::  vpt_surface             !< near-surface virtual potential temperature
718       REAL(wp) ::  zi_l                    !< mean boundary-layer depth on subdomain   
719       REAL(wp) ::  zi_local                !< local boundary-layer depth 
720   
721       REAL(wp), DIMENSION(nzb:nzt+1) ::  vpt_col !< vertical profile of virtual potential temperature at (j,i)-grid point
722
723       
724!
725!--    Calculate mean boundary-layer height from boundary data.
726!--    Start with the left and right boundaries.
727       zi_l      = 0.0_wp
728       IF ( bc_dirichlet_l  .OR.  bc_dirichlet_r )  THEN
729!
730!--       Determine index along x. Please note, index indicates boundary
731!--       grid point for scalars.
732          i = MERGE( -1, nxr + 1, bc_dirichlet_l )
733   
734          DO  j = nys, nyn
735!
736!--          Determine topography top index at current (j,i) index
737             k_surface = get_topography_top_index_ji( j, i, 's' )
738!
739!--          Pre-compute surface virtual temperature. Therefore, use 2nd
740!--          prognostic level according to Heinze et al. (2017).
741             IF ( humidity )  THEN
742                vpt_surface = pt(k_surface+2,j,i) *                            &
743                            ( 1.0_wp + 0.61_wp * q(k_surface+2,j,i) )
744                vpt_col     = pt(:,j,i) * ( 1.0_wp + 0.61_wp * q(:,j,i) )
745             ELSE
746                vpt_surface = pt(k_surface+2,j,i)
747                vpt_col     = pt(:,j,i)
748             ENDIF
749!
750!--          Calculate local boundary layer height from bulk Richardson number,
751!--          i.e. the height where the bulk Richardson number exceeds its
752!--          critical value of 0.25 (according to Heinze et al., 2017).
753!--          Note, no interpolation of u- and v-component is made, as both
754!--          are mainly mean inflow profiles with very small spatial variation.
755             zi_local = 0.0_wp
756             DO  k = k_surface+1, nzt
757                u_comp = MERGE( u(k,j,i+1), u(k,j,i), bc_dirichlet_l )
758                v_comp = v(k,j,i)
759                ri_bulk = zu(k) * g / vpt_surface *                            &
760                          ( vpt_col(k) - vpt_surface ) /                       &
761                          ( u_comp * u_comp + v_comp * v_comp ) 
762                       
763                IF ( zi_local == 0.0_wp  .AND.  ri_bulk > ri_bulk_crit )       &
764                   zi_local = zu(k)           
765             ENDDO
766!
767!--          Assure that the minimum local boundary-layer depth is at least at
768!--          the second vertical grid level.
769             zi_l = zi_l + MAX( zi_local, zu(k_surface+2) )   
770             
771          ENDDO
772       
773       ENDIF
774!
775!--    Do the same at the north and south boundaries.
776       IF ( bc_dirichlet_s  .OR.  bc_dirichlet_n )  THEN
777       
778          j = MERGE( -1, nyn + 1, bc_dirichlet_s )
779       
780          DO  i = nxl, nxr
781             k_surface = get_topography_top_index_ji( j, i, 's' )
782 
783             IF ( humidity )  THEN
784                vpt_surface = pt(k_surface+2,j,i) *                            &
785                            ( 1.0_wp + 0.61_wp * q(k_surface+2,j,i) )
786                vpt_col     = pt(:,j,i) * ( 1.0_wp + 0.61_wp * q(:,j,i) )
787             ELSE
788                vpt_surface = pt(k_surface+2,j,i)
789                vpt_col  = pt(:,j,i)
790             ENDIF
791         
792             zi_local = 0.0_wp
793             DO  k = k_surface+1, nzt               
794                u_comp = u(k,j,i)
795                v_comp = MERGE( v(k,j+1,i), v(k,j,i), bc_dirichlet_s )
796                ri_bulk = zu(k) * g / vpt_surface *                            &
797                       ( vpt_col(k) - vpt_surface ) /                          &
798                       ( u_comp * u_comp + v_comp * v_comp ) 
799                       
800                IF ( zi_local == 0.0_wp  .AND.  ri_bulk > 0.25_wp )            &
801                   zi_local = zu(k)   
802             ENDDO
803             zi_l = zi_l + MAX( zi_local, zu(k_surface+2) )   
804         
805          ENDDO
806         
807       ENDIF
808   
809#if defined( __parallel )
810       CALL MPI_ALLREDUCE( zi_l, zi_ribulk, 1, MPI_REAL, MPI_SUM,              &
811                           comm2d, ierr )
812#else
813       zi_ribulk = zi_l
814#endif
815       zi_ribulk = zi_ribulk / REAL( 2 * nx + 2 * ny, KIND = wp )
816!
817!--    Finally, check if boundary layer depth is not below the any topography.
818!--    zi_ribulk will be used to adjust rayleigh damping height, i.e. the
819!--    lower level of the sponge layer, as well as to adjust the synthetic
820!--    turbulence generator accordingly. If Rayleigh damping would be applied
821!--    near buildings, etc., this would spoil the simulation results.
822       topo_max_l = zw(MAXVAL( get_topography_top_index( 's' )))
823       
824#if defined( __parallel )
825       CALL MPI_ALLREDUCE( topo_max_l, topo_max, 1, MPI_REAL, MPI_MAX,            &
826                           comm2d, ierr )
827#else
828       topo_max     = topo_max_l
829#endif
830
831       zi_ribulk = MAX( zi_ribulk, topo_max )
832
833    END SUBROUTINE calc_zi
834   
835   
836!------------------------------------------------------------------------------!
837! Description:
838!------------------------------------------------------------------------------!
839!> Adjust the height where the rayleigh damping starts, i.e. the lower level
840!> of the sponge layer.
841!------------------------------------------------------------------------------!
842    SUBROUTINE adjust_sponge_layer
843       
844       USE arrays_3d,                                                          &
845           ONLY:  rdf, rdf_sc, zu
846       
847       USE basic_constants_and_equations_mod,                                  &
848           ONLY:  pi
849           
850       USE control_parameters,                                                 &
851           ONLY:  rayleigh_damping_height, rayleigh_damping_factor
852       
853       USE kinds
854
855       IMPLICIT NONE
856
857       INTEGER(iwp) :: k   !< loop index in z-direction
858   
859       REAL(wp) ::  rdh    !< updated Rayleigh damping height
860 
861   
862       IF ( rayleigh_damping_height > 0.0_wp  .AND.                            &
863            rayleigh_damping_factor > 0.0_wp )  THEN
864!
865!--       Update Rayleigh-damping height and re-calculate height-depending
866!--       damping coefficients.
867!--       Assure that rayleigh damping starts well above the boundary layer.   
868          rdh = MIN( MAX( zi_ribulk * 1.3_wp, 10.0_wp * dz(1) ),               & 
869                     0.8_wp * zu(nzt), rayleigh_damping_height )
870!
871!--       Update Rayleigh damping factor
872          DO  k = nzb+1, nzt
873             IF ( zu(k) >= rdh )  THEN
874                rdf(k) = rayleigh_damping_factor *                             &
875                                          ( SIN( pi * 0.5_wp * ( zu(k) - rdh ) &
876                                        / ( zu(nzt) - rdh ) )                  &
877                                          )**2
878             ENDIF
879          ENDDO
880          rdf_sc = rdf
881       
882       ENDIF
883
884    END SUBROUTINE adjust_sponge_layer
885   
886!------------------------------------------------------------------------------!
887! Description:
888! ------------
889!> Performs consistency checks
890!------------------------------------------------------------------------------!
891    SUBROUTINE nesting_offl_check_parameters
892   
893       USE control_parameters,                                                    &
894        ONLY:  child_domain, message_string, nesting_offline
895
896       IMPLICIT NONE
897!
898!--    Perform checks
899       IF ( nesting_offline  .AND.  child_domain )  THEN
900          message_string = 'Offline nesting is only applicable in root model.'
901          CALL message( 'stg_check_parameters', 'PA0622', 1, 2, 0, 6, 0 )       
902       ENDIF
903
904
905    END SUBROUTINE nesting_offl_check_parameters
906   
907!------------------------------------------------------------------------------!
908! Description:
909! ------------
910!> Reads the parameter list nesting_offl_parameters
911!------------------------------------------------------------------------------!
912    SUBROUTINE nesting_offl_parin
913
914       IMPLICIT NONE
915       
916       CHARACTER (LEN=80) ::  line   !< dummy string that contains the current line of the parameter file
917
918
919       NAMELIST /nesting_offl_parameters/   nesting_offline
920
921       line = ' '
922
923!
924!--    Try to find stg package
925       REWIND ( 11 )
926       line = ' '
927       DO WHILE ( INDEX( line, '&nesting_offl_parameters' ) == 0 )
928          READ ( 11, '(A)', END=20 )  line
929       ENDDO
930       BACKSPACE ( 11 )
931
932!
933!--    Read namelist
934       READ ( 11, nesting_offl_parameters, ERR = 10, END = 20 )
935
936       GOTO 20
937
938 10    BACKSPACE( 11 )
939       READ( 11 , '(A)') line
940       CALL parin_fail_message( 'nesting_offl_parameters', line )
941
942 20    CONTINUE
943
944
945    END SUBROUTINE nesting_offl_parin
946
947!------------------------------------------------------------------------------!
948! Description:
949! ------------
950!> Writes information about offline nesting into HEADER file
951!------------------------------------------------------------------------------!
952    SUBROUTINE nesting_offl_header ( io )
953
954       IMPLICIT NONE
955
956       INTEGER(iwp), INTENT(IN) ::  io !< Unit of the output file
957
958       WRITE ( io, 1 )
959       IF ( nesting_offline )  THEN
960          WRITE ( io, 3 )
961       ELSE
962          WRITE ( io, 2 )
963       ENDIF
964
9651 FORMAT (//' Offline nesting in COSMO model:'/                                &
966              ' -------------------------------'/)
9672 FORMAT (' --> No offlince nesting is used (default) ')
9683 FORMAT (' --> Offlince nesting is used. Boundary data is read from dynamic input file ')
969
970    END SUBROUTINE nesting_offl_header
971
972!------------------------------------------------------------------------------!
973! Description:
974! ------------
975!> Allocate arrays used to read boundary data from NetCDF file and initialize
976!> boundary data.
977!------------------------------------------------------------------------------!
978    SUBROUTINE nesting_offl_init
979   
980       USE netcdf_data_input_mod,                                              &
981           ONLY:  netcdf_data_input_offline_nesting 
982
983       IMPLICIT NONE
984
985
986!--    Allocate arrays for geostrophic wind components. Arrays will
987!--    incorporate 2 time levels in order to interpolate in between.
988       ALLOCATE( nest_offl%ug(0:1,1:nzt) )
989       ALLOCATE( nest_offl%vg(0:1,1:nzt) )
990!
991!--    Allocate arrays for reading boundary values. Arrays will incorporate 2
992!--    time levels in order to interpolate in between.
993       IF ( bc_dirichlet_l )  THEN
994          ALLOCATE( nest_offl%u_left(0:1,nzb+1:nzt,nys:nyn)  )
995          ALLOCATE( nest_offl%v_left(0:1,nzb+1:nzt,nysv:nyn) )
996          ALLOCATE( nest_offl%w_left(0:1,nzb+1:nzt-1,nys:nyn) )
997          IF ( humidity )       ALLOCATE( nest_offl%q_left(0:1,nzb+1:nzt,nys:nyn)  )
998          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_left(0:1,nzb+1:nzt,nys:nyn) )
999       ENDIF
1000       IF ( bc_dirichlet_r )  THEN
1001          ALLOCATE( nest_offl%u_right(0:1,nzb+1:nzt,nys:nyn)  )
1002          ALLOCATE( nest_offl%v_right(0:1,nzb+1:nzt,nysv:nyn) )
1003          ALLOCATE( nest_offl%w_right(0:1,nzb+1:nzt-1,nys:nyn) )
1004          IF ( humidity )       ALLOCATE( nest_offl%q_right(0:1,nzb+1:nzt,nys:nyn)  )
1005          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_right(0:1,nzb+1:nzt,nys:nyn) )
1006       ENDIF
1007       IF ( bc_dirichlet_n )  THEN
1008          ALLOCATE( nest_offl%u_north(0:1,nzb+1:nzt,nxlu:nxr) )
1009          ALLOCATE( nest_offl%v_north(0:1,nzb+1:nzt,nxl:nxr)  )
1010          ALLOCATE( nest_offl%w_north(0:1,nzb+1:nzt-1,nxl:nxr) )
1011          IF ( humidity )       ALLOCATE( nest_offl%q_north(0:1,nzb+1:nzt,nxl:nxr)  )
1012          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_north(0:1,nzb+1:nzt,nxl:nxr) )
1013       ENDIF
1014       IF ( bc_dirichlet_s )  THEN
1015          ALLOCATE( nest_offl%u_south(0:1,nzb+1:nzt,nxlu:nxr) )
1016          ALLOCATE( nest_offl%v_south(0:1,nzb+1:nzt,nxl:nxr)  )
1017          ALLOCATE( nest_offl%w_south(0:1,nzb+1:nzt-1,nxl:nxr)    )
1018          IF ( humidity )       ALLOCATE( nest_offl%q_south(0:1,nzb+1:nzt,nxl:nxr)  )
1019          IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_south(0:1,nzb+1:nzt,nxl:nxr) )
1020       ENDIF
1021       
1022       ALLOCATE( nest_offl%u_top(0:1,nys:nyn,nxlu:nxr) )
1023       ALLOCATE( nest_offl%v_top(0:1,nysv:nyn,nxl:nxr) )
1024       ALLOCATE( nest_offl%w_top(0:1,nys:nyn,nxl:nxr)  )
1025       IF ( humidity )       ALLOCATE( nest_offl%q_top(0:1,nys:nyn,nxl:nxr)  )
1026       IF ( .NOT. neutral )  ALLOCATE( nest_offl%pt_top(0:1,nys:nyn,nxl:nxr) )
1027
1028!
1029!--    Read COSMO data at lateral and top boundaries
1030       CALL netcdf_data_input_offline_nesting
1031!
1032!--    Initialize boundary data.
1033       IF ( bc_dirichlet_l )  THEN
1034          u(nzb+1:nzt,nys:nyn,0)    = nest_offl%u_left(0,nzb+1:nzt,nys:nyn)
1035          v(nzb+1:nzt,nysv:nyn,-1)  = nest_offl%v_left(0,nzb+1:nzt,nysv:nyn) 
1036          w(nzb+1:nzt-1,nys:nyn,-1) = nest_offl%w_left(0,nzb+1:nzt-1,nys:nyn)
1037          IF ( .NOT. neutral )  pt(nzb+1:nzt,nys:nyn,-1) =                     &
1038                                   nest_offl%pt_left(0,nzb+1:nzt,nys:nyn)
1039          IF ( humidity      )  q(nzb+1:nzt,nys:nyn,-1)  =                     &
1040                                   nest_offl%q_left(0,nzb+1:nzt,nys:nyn) 
1041       ENDIF
1042       IF ( bc_dirichlet_r )  THEN
1043          u(nzb+1:nzt,nys:nyn,nxr+1)   = nest_offl%u_right(0,nzb+1:nzt,nys:nyn) 
1044          v(nzb+1:nzt,nysv:nyn,nxr+1)  = nest_offl%v_right(0,nzb+1:nzt,nysv:nyn)
1045          w(nzb+1:nzt-1,nys:nyn,nxr+1) = nest_offl%w_right(0,nzb+1:nzt-1,nys:nyn)
1046          IF ( .NOT. neutral )  pt(nzb+1:nzt,nys:nyn,nxr+1) =                  &
1047                                   nest_offl%pt_right(0,nzb+1:nzt,nys:nyn)
1048          IF ( humidity      )  q(nzb+1:nzt,nys:nyn,nxr+1)  =                  &
1049                                   nest_offl%q_right(0,nzb+1:nzt,nys:nyn)
1050       ENDIF
1051       IF ( bc_dirichlet_s )  THEN
1052          u(nzb+1:nzt,-1,nxlu:nxr)  = nest_offl%u_south(0,nzb+1:nzt,nxlu:nxr)
1053          v(nzb+1:nzt,0,nxl:nxr)    = nest_offl%v_south(0,nzb+1:nzt,nxl:nxr) 
1054          w(nzb+1:nzt-1,-1,nxl:nxr) = nest_offl%w_south(0,nzb+1:nzt-1,nxl:nxr) 
1055          IF ( .NOT. neutral )  pt(nzb+1:nzt,-1,nxl:nxr) =                     &
1056                                   nest_offl%pt_south(0,nzb+1:nzt,nxl:nxr)
1057          IF ( humidity      )  q(nzb+1:nzt,-1,nxl:nxr)  =                     &
1058                                   nest_offl%q_south(0,nzb+1:nzt,nxl:nxr) 
1059       ENDIF
1060       IF ( bc_dirichlet_n )  THEN
1061          u(nzb+1:nzt,nyn+1,nxlu:nxr)  = nest_offl%u_north(0,nzb+1:nzt,nxlu:nxr)
1062          v(nzb+1:nzt,nyn+1,nxl:nxr)   = nest_offl%v_north(0,nzb+1:nzt,nxl:nxr) 
1063          w(nzb+1:nzt-1,nyn+1,nxl:nxr) = nest_offl%w_north(0,nzb+1:nzt-1,nxl:nxr)
1064          IF ( .NOT. neutral )  pt(nzb+1:nzt,nyn+1,nxl:nxr) =                  &
1065                                   nest_offl%pt_north(0,nzb+1:nzt,nxl:nxr) 
1066          IF ( humidity      )  q(nzb+1:nzt,nyn+1,nxl:nxr)  =                  &
1067                                   nest_offl%q_north(0,nzb+1:nzt,nxl:nxr)
1068       ENDIF
1069!
1070!--    Initialize geostrophic wind components. Actually this is already done in
1071!--    init_3d_model when initializing_action = 'inifor', however, in speical
1072!--    case of user-defined initialization this will be done here again, in
1073!--    order to have a consistent initialization.
1074       ug(nzb+1:nzt) = nest_offl%ug(0,nzb+1:nzt)
1075       vg(nzb+1:nzt) = nest_offl%vg(0,nzb+1:nzt)
1076!
1077!--    Set bottom and top boundary condition for geostrophic wind components
1078       ug(nzt+1) = ug(nzt)
1079       vg(nzt+1) = vg(nzt)
1080       ug(nzb)   = ug(nzb+1)
1081       vg(nzb)   = vg(nzb+1)
1082!
1083!--    Initial calculation of the boundary layer depth from the prescribed
1084!--    boundary data. This is requiered for initialize the synthetic turbulence
1085!--    generator correctly.
1086       CALL calc_zi
1087       
1088!
1089!--    After boundary data is initialized, mask topography at the
1090!--    boundaries for the velocity components.
1091       u = MERGE( u, 0.0_wp, BTEST( wall_flags_0, 1 ) )
1092       v = MERGE( v, 0.0_wp, BTEST( wall_flags_0, 2 ) )
1093       w = MERGE( w, 0.0_wp, BTEST( wall_flags_0, 3 ) )
1094       
1095       CALL calc_zi
1096       
1097!
1098!--    After boundary data is initialized, ensure mass conservation
1099       CALL nesting_offl_mass_conservation
1100
1101    END SUBROUTINE nesting_offl_init
1102   
1103!------------------------------------------------------------------------------!
1104! Description:
1105!------------------------------------------------------------------------------!
1106!> Interpolation function, used to interpolate boundary data in time.
1107!------------------------------------------------------------------------------!
1108    FUNCTION interpolate_in_time( var_t1, var_t2, fac  ) 
1109       
1110       USE kinds
1111
1112       IMPLICIT NONE
1113
1114       REAL(wp)            :: interpolate_in_time !< time-interpolated boundary value
1115       REAL(wp)            :: var_t1              !< boundary value at t1
1116       REAL(wp)            :: var_t2              !< boundary value at t2
1117       REAL(wp)            :: fac                 !< interpolation factor
1118
1119       interpolate_in_time = ( 1.0_wp - fac ) * var_t1 + fac * var_t2     
1120
1121    END FUNCTION interpolate_in_time
1122
1123
1124
1125 END MODULE nesting_offl_mod
Note: See TracBrowser for help on using the repository browser.