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

Last change on this file since 3736 was 3705, checked in by suehring, 6 years ago

last commit documented

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