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

Last change on this file since 3885 was 3885, checked in by kanani, 2 years ago

restructure/add location/debug messages

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