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

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

Consider time-dependent geostrophic wind components in offline nesting

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