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

Last change on this file since 3382 was 3364, checked in by suehring, 6 years ago

Set svn id for nesting_offl_mod

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