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

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

svn keyword forgotten to set for nesting_offl_mod

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