source: palm/trunk/SOURCE/large_scale_forcing_nudging_mod.f90 @ 3183

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

last commit documented

  • Property svn:keywords set to Id
File size: 78.2 KB
Line 
1!> @file large_scale_forcing_nudging_mod.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2018 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: large_scale_forcing_nudging_mod.f90 3183 2018-07-27 14:25:55Z suehring $
27! * Adjustment to new Inifor version:
28!   - No vertical interpolation/extrapolation of lateral boundary data required
29!     any more (Inifor can treat grid stretching now
30!   - Revise initialization in case of COSMO forcing
31! * Rename variables and subroutines for offline nesting
32!
33! 3182 2018-07-27 13:36:03Z suehring
34! Error messages revised
35!
36! 3045 2018-05-28 07:55:41Z Giersch
37! Error messages revised
38!
39! 3026 2018-05-22 10:30:53Z schwenkel
40! Changed the name specific humidity to mixing ratio, since we are computing
41! mixing ratios.
42!
43! 2970 2018-04-13 15:09:23Z suehring
44! Bugfix in old large-scale forcing mode
45!
46! 2938 2018-03-27 15:52:42Z suehring
47! Further improvements for nesting in larger-scale model
48!
49! 2863 2018-03-08 11:36:25Z suehring
50! Corrected "Former revisions" section
51!
52! 2696 2017-12-14 17:12:51Z kanani
53! Change in file header (GPL part)
54! Forcing with larger-scale models implemented (MS)
55!
56! 2342 2017-08-08 11:00:43Z boeske
57! fixed check if surface forcing data is available until end of simulation
58!
59! 2320 2017-07-21 12:47:43Z suehring
60! initial revision
61!
62! Description:
63! ------------
64!> Calculates large scale forcings (geostrophic wind and subsidence velocity) as
65!> well as surfaces fluxes dependent on time given in an external file (LSF_DATA).
66!> Moreover, module contains nudging routines, where u, v, pt and q are nudged
67!> to given profiles on a relaxation timescale tnudge.
68!> Profiles are read in from NUDGING_DATA.
69!> Code is based on Neggers et al. (2012) and also in parts on DALES and UCLA-LES.
70!> @todo: Revise reading of ASCII-files
71!> @todo: Remove unused variables and control flags
72!> @todo: Revise large-scale facing of surface variables
73!> @todo: Revise control flags lsf_exception, lsf_surf, lsf_vert, etc.
74!--------------------------------------------------------------------------------!
75 MODULE lsf_nudging_mod
76
77    USE arrays_3d,                                                             &
78        ONLY:  dzw, e, diss, heatflux_input_conversion, pt, pt_init, q,        &
79               q_init, s, tend, u, u_init, ug, v, v_init, vg, w, w_subs,       &
80               waterflux_input_conversion, zu, zw                 
81
82    USE control_parameters,                                                    &
83        ONLY:  bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, bc_dirichlet_s, &
84               bc_lr, bc_ns, bc_pt_b, bc_q_b, constant_diffusion,              &
85               constant_heatflux, constant_waterflux,                          &
86               data_output_pr, dt_3d, end_time,                                &
87               humidity, initializing_actions, intermediate_timestep_count,    &
88               ibc_pt_b, ibc_q_b,                                              &
89               large_scale_forcing, large_scale_subsidence, lsf_surf, lsf_vert,&
90               lsf_exception, message_string, nesting_offline, neutral,        &
91               nudging, passive_scalar, pt_surface, ocean, q_surface,          &
92               surface_heatflux, surface_pressure, surface_waterflux,          &
93               topography, use_subsidence_tendencies
94
95    USE grid_variables
96
97    USE indices,                                                               &
98        ONLY:  nbgp, ngp_sums_ls, nx, nxl, nxlg, nxlu, nxr, nxrg, ny, nys,     &
99               nysv, nysg, nyn, nyng, nzb, nz, nzt, wall_flags_0
100
101    USE kinds
102
103    USE netcdf_data_input_mod,                                                 &
104        ONLY:  nest_offl
105
106    USE pegrid
107
108    USE surface_mod,                                                           &
109        ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
110
111    USE statistics,                                                            &
112        ONLY:  hom, statistic_regions, sums_ls_l, weight_substep
113
114    INTEGER(iwp) ::  nlsf = 1000                       !< maximum number of profiles in LSF_DATA (large scale forcing)
115    INTEGER(iwp) ::  ntnudge = 1000                    !< maximum number of profiles in NUDGING_DATA (nudging)
116
117    REAL(wp) ::  d_area_t
118
119    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ptnudge     !< vertical profile of pot. temperature interpolated to vertical grid (nudging)
120    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  qnudge      !< vertical profile of water vapor mixing ratio interpolated to vertical grid (nudging)
121    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tnudge      !< vertical profile of nudging time scale interpolated to vertical grid (nudging) 
122    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  td_lsa_lpt  !< temperature tendency due to large scale advection (large scale forcing)
123    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  td_lsa_q    !< water vapor mixing ratio tendency due to large scale advection (large scale forcing)
124    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  td_sub_lpt  !< temperature tendency due to subsidence/ascent (large scale forcing)
125    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  td_sub_q    !< water vapor mixing ratio tendency due to subsidence/ascent (large scale forcing)
126    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ug_vert     !< vertical profile of geostrophic wind component in x-direction interpolated to vertical grid (large scale forcing)
127    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  unudge      !< vertical profile of wind component in x-direction interpolated to vertical grid (nudging)
128    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  vnudge      !< vertical profile of wind component in y-direction interpolated to vertical grid (nudging)
129    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  vg_vert     !< vertical profile of geostrophic wind component in y-direction interpolated to vertical grid (large scale forcing)
130    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  wnudge      !< vertical profile of subsidence/ascent velocity interpolated to vertical grid (nudging) ???
131    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  wsubs_vert  !< vertical profile of wind component in z-direction interpolated to vertical grid (nudging) ???
132
133    REAL(wp), DIMENSION(:), ALLOCATABLE ::  shf_surf      !< time-dependent surface sensible heat flux (large scale forcing)
134    REAL(wp), DIMENSION(:), ALLOCATABLE ::  timenudge     !< times at which vertical profiles are defined in NUDGING_DATA (nudging)
135    REAL(wp), DIMENSION(:), ALLOCATABLE ::  time_surf     !< times at which surface values/fluxes are defined in LSF_DATA (large scale forcing)
136    REAL(wp), DIMENSION(:), ALLOCATABLE ::  time_vert     !< times at which vertical profiles are defined in LSF_DATA (large scale forcing)
137
138    REAL(wp), DIMENSION(:), ALLOCATABLE ::  tmp_tnudge    !< current nudging time scale
139
140    REAL(wp), DIMENSION(:), ALLOCATABLE ::  p_surf        !< time-dependent surface pressure (large scale forcing)
141    REAL(wp), DIMENSION(:), ALLOCATABLE ::  pt_surf       !< time-dependent surface temperature (large scale forcing)
142    REAL(wp), DIMENSION(:), ALLOCATABLE ::  qsws_surf     !< time-dependent surface latent heat flux (large scale forcing)
143    REAL(wp), DIMENSION(:), ALLOCATABLE ::  q_surf        !< time-dependent surface water vapor mixing ratio (large scale forcing)
144
145    SAVE
146    PRIVATE
147!
148!-- Public subroutines
149    PUBLIC calc_tnudge, ls_forcing_surf, ls_forcing_vert, ls_advec, lsf_init,  &
150           lsf_nudging_check_parameters, nudge_init,                           &
151           lsf_nudging_check_data_output_pr, lsf_nudging_header,               &
152           lsf_nesting_offline, lsf_nesting_offline_mass_conservation,         & 
153           nudge, nudge_ref
154           
155!
156!-- Public variables
157    PUBLIC qsws_surf, shf_surf, td_lsa_lpt, td_lsa_q, td_sub_lpt,              &
158           td_sub_q, time_vert
159
160
161    INTERFACE ls_advec
162       MODULE PROCEDURE ls_advec
163       MODULE PROCEDURE ls_advec_ij
164    END INTERFACE ls_advec
165
166    INTERFACE nudge
167       MODULE PROCEDURE nudge
168       MODULE PROCEDURE nudge_ij
169    END INTERFACE nudge
170
171 CONTAINS
172
173
174!------------------------------------------------------------------------------!
175! Description:
176! ------------
177!> In this subroutine a constant mass within the model domain is guaranteed.
178!> Larger-scale models may be based on a compressible equation system, which is
179!> not consistent with PALMs incompressible equation system. In order to avoid
180!> a decrease or increase of mass during the simulation, non-divergent flow
181!> through the lateral and top boundaries is compensated by the vertical wind
182!> component at the top boundary.
183!------------------------------------------------------------------------------!
184    SUBROUTINE lsf_nesting_offline_mass_conservation
185
186       USE control_parameters,                                                 &
187           ONLY:  volume_flow
188
189       IMPLICIT NONE
190
191       INTEGER(iwp) ::  i !< grid index in x-direction
192       INTEGER(iwp) ::  j !< grid index in y-direction
193       INTEGER(iwp) ::  k !< grid index in z-direction
194
195       REAL(wp) ::  w_correct                       !< vertical velocity increment required to compensate non-divergent flow through the boundaries
196       REAL(wp), DIMENSION(1:3) ::  volume_flow_l   !< local volume flow
197
198       volume_flow   = 0.0_wp
199       volume_flow_l = 0.0_wp
200
201       d_area_t = 1.0_wp / ( ( nx + 1 ) * dx * ( ny + 1 ) * dy )
202
203       IF ( bc_dirichlet_l )  THEN
204          i = nxl
205          DO  j = nys, nyn
206             DO  k = nzb+1, nzt
207                volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k) * dy   &
208                                     * MERGE( 1.0_wp, 0.0_wp,                  &
209                                              BTEST( wall_flags_0(k,j,i), 1 ) )
210             ENDDO
211          ENDDO
212       ENDIF
213       IF ( bc_dirichlet_r )  THEN
214          i = nxr+1
215          DO  j = nys, nyn
216             DO  k = nzb+1, nzt
217                volume_flow_l(1) = volume_flow_l(1) - u(k,j,i) * dzw(k) * dy   &
218                                     * MERGE( 1.0_wp, 0.0_wp,                  &
219                                              BTEST( wall_flags_0(k,j,i), 1 ) )
220             ENDDO
221          ENDDO
222       ENDIF
223       IF ( bc_dirichlet_s )  THEN
224          j = nys
225          DO  i = nxl, nxr
226             DO  k = nzb+1, nzt
227                volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k) * dx   &
228                                     * MERGE( 1.0_wp, 0.0_wp,                  &
229                                              BTEST( wall_flags_0(k,j,i), 2 ) )
230             ENDDO
231          ENDDO
232       ENDIF
233       IF ( bc_dirichlet_n )  THEN
234          j = nyn+1
235          DO  i = nxl, nxr
236             DO  k = nzb+1, nzt
237                volume_flow_l(2) = volume_flow_l(2) - v(k,j,i) * dzw(k) * dx   &
238                                     * MERGE( 1.0_wp, 0.0_wp,                  &
239                                              BTEST( wall_flags_0(k,j,i), 2 ) )
240             ENDDO
241          ENDDO
242       ENDIF
243!
244!--    Top boundary
245       k = nzt
246       DO  i = nxl, nxr
247          DO  j = nys, nyn
248             volume_flow_l(3) = volume_flow_l(3) - w(k,j,i) * dx * dy
249          ENDDO
250       ENDDO
251
252#if defined( __parallel )
253       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
254       CALL MPI_ALLREDUCE( volume_flow_l, volume_flow, 3, MPI_REAL, MPI_SUM,      &
255                           comm2d, ierr )
256#else
257       volume_flow = volume_flow_l
258#endif
259
260       w_correct = SUM( volume_flow ) * d_area_t
261
262       DO  i = nxl, nxr
263          DO  j = nys, nyn
264             DO  k = nzt, nzt + 1
265                w(k,j,i) = w(k,j,i) + w_correct
266             ENDDO
267          ENDDO
268       ENDDO
269
270    END SUBROUTINE lsf_nesting_offline_mass_conservation
271
272
273!------------------------------------------------------------------------------!
274! Description:
275! ------------
276!> Set the lateral and top boundary conditions in case the PALM domain is
277!> nested offline in a mesoscale model.
278!------------------------------------------------------------------------------!
279    SUBROUTINE lsf_nesting_offline
280
281       USE control_parameters,                                                 &
282           ONLY:  bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r,              &
283                  bc_dirichlet_s, humidity, neutral, passive_scalar, rans_mode,&
284                  rans_tke_e, time_since_reference_point                     
285
286       IMPLICIT NONE
287
288       INTEGER(iwp) ::  i !< running index x-direction
289       INTEGER(iwp) ::  j !< running index y-direction
290       INTEGER(iwp) ::  k !< running index z-direction
291       INTEGER(iwp) ::  t !< running index for time levels
292
293       REAL(wp) ::  ddt_lsf !< inverse value of time resolution of forcing data
294       REAL(wp) ::  t_ref   !< time past since last reference step
295     
296!
297!--    Calculate time interval of forcing data       
298       ddt_lsf = 1.0_wp / ( nest_offl%time(nest_offl%tind_p) -                 &
299                            nest_offl%time(nest_offl%tind) )
300!
301!--    Calculate reziproke time past since last reference step. Please note,
302!--    the time coordinate is still not updated, so that the actual time need
303!--    to be incremented by dt_3d. Moreover, note that the simulation time
304!--    passed since simulation start is time_since_reference_point, not
305!--    simulated_time!
306       t_ref = time_since_reference_point + dt_3d -                            &
307                                            nest_offl%time(nest_offl%tind)
308                                           
309       IF ( bc_dirichlet_l )  THEN
310
311          DO  j = nys, nyn
312             DO  k = nzb+1, nzt
313                u(k,j,nxlg:nxl) = nest_offl%u_left(0,k,j) + ddt_lsf * t_ref  * &
314                       ( nest_offl%u_left(1,k,j) - nest_offl%u_left(0,k,j) ) * &
315                         MERGE( 1.0_wp, 0.0_wp,                                &
316                                BTEST( wall_flags_0(k,j,nxlg:nxl), 1 ) )
317             ENDDO
318          ENDDO
319
320          DO  j = nys, nyn
321             DO  k = nzb+1, nzt-1
322                w(k,j,nxlg:nxl-1) = nest_offl%w_left(0,k,j) + ddt_lsf * t_ref *&
323                       ( nest_offl%w_left(1,k,j) - nest_offl%w_left(0,k,j) )  *&
324                         MERGE( 1.0_wp, 0.0_wp,                                &
325                                BTEST( wall_flags_0(k,j,nxlg:nxl-1), 3 ) )
326             ENDDO
327          ENDDO
328
329          DO  j = nysv, nyn
330             DO  k = nzb+1, nzt
331                v(k,j,nxlg:nxl-1) = nest_offl%v_left(0,k,j) + ddt_lsf * t_ref *&
332                       ( nest_offl%v_left(1,k,j) - nest_offl%v_left(0,k,j) )  *&
333                         MERGE( 1.0_wp, 0.0_wp,                                &
334                                BTEST( wall_flags_0(k,j,nxlg:nxl-1), 2 ) )
335             ENDDO
336          ENDDO
337
338          IF ( .NOT. neutral )  THEN
339             DO  j = nys, nyn
340                DO  k = nzb+1, nzt
341                   pt(k,j,nxlg:nxl-1) = nest_offl%pt_left(0,k,j) + ddt_lsf *   &
342                                                                   t_ref   *   &
343                       ( nest_offl%pt_left(1,k,j) - nest_offl%pt_left(0,k,j) )
344 
345                ENDDO
346             ENDDO
347          ENDIF
348
349          IF ( humidity )  THEN
350             DO  j = nys, nyn
351                DO  k = nzb+1, nzt
352                   q(k,j,nxlg:nxl-1) = nest_offl%q_left(0,k,j) + ddt_lsf   *   &
353                                                                 t_ref     *   &
354                       ( nest_offl%q_left(1,k,j) - nest_offl%q_left(0,k,j) )
355 
356                ENDDO
357             ENDDO
358          ENDIF
359
360       ENDIF
361
362       IF ( bc_dirichlet_r )  THEN
363
364          DO  j = nys, nyn
365             DO  k = nzb+1, nzt
366                u(k,j,nxr+1:nxrg) = nest_offl%u_right(0,k,j) + ddt_lsf * t_ref *&
367                      ( nest_offl%u_right(1,k,j) - nest_offl%u_right(0,k,j) )  *&
368                         MERGE( 1.0_wp, 0.0_wp,                                 &
369                                BTEST( wall_flags_0(k,j,nxr+1:nxrg), 1 ) )
370
371             ENDDO
372          ENDDO
373          DO  j = nys, nyn
374             DO  k = nzb+1, nzt-1
375                w(k,j,nxr+1:nxrg) = nest_offl%w_right(0,k,j) + ddt_lsf * t_ref *&
376                      ( nest_offl%w_right(1,k,j) - nest_offl%w_right(0,k,j) )  *&
377                         MERGE( 1.0_wp, 0.0_wp,                                 &
378                                BTEST( wall_flags_0(k,j,nxr+1:nxrg), 3 ) )
379             ENDDO
380          ENDDO
381
382          DO  j = nysv, nyn
383             DO  k = nzb+1, nzt
384                v(k,j,nxr+1:nxrg) = nest_offl%v_right(0,k,j) + ddt_lsf * t_ref *&
385                      ( nest_offl%v_right(1,k,j) - nest_offl%v_right(0,k,j) )  *&
386                         MERGE( 1.0_wp, 0.0_wp,                                 &
387                                BTEST( wall_flags_0(k,j,nxr+1:nxrg), 2 ) )
388             ENDDO
389          ENDDO
390
391          IF ( .NOT. neutral )  THEN
392             DO  j = nys, nyn
393                DO  k = nzb+1, nzt
394                   pt(k,j,nxr+1:nxrg) = nest_offl%pt_right(0,k,j) + ddt_lsf *  &
395                                                                    t_ref   *  &
396                     ( nest_offl%pt_right(1,k,j) - nest_offl%pt_right(0,k,j) )
397 
398                ENDDO
399             ENDDO
400          ENDIF
401
402          IF ( humidity )  THEN
403             DO  j = nys, nyn
404                DO  k = nzb+1, nzt
405                   q(k,j,nxr+1:nxrg) = nest_offl%q_right(0,k,j) + ddt_lsf   *  &
406                                                                    t_ref   *  &
407                       ( nest_offl%q_right(1,k,j) - nest_offl%q_right(0,k,j) )
408 
409                ENDDO
410             ENDDO
411          ENDIF
412
413       ENDIF
414
415       IF ( bc_dirichlet_s )  THEN
416
417          DO  i = nxl, nxr
418             DO  k = nzb+1, nzt
419                v(k,nysg:nys,i)   = nest_offl%v_south(0,k,i) + ddt_lsf * t_ref *&
420                      ( nest_offl%v_south(1,k,i) - nest_offl%v_south(0,k,i) )  *&
421                         MERGE( 1.0_wp, 0.0_wp,                                 &
422                                BTEST( wall_flags_0(k,nysg:nys,i), 2 ) )
423             ENDDO
424          ENDDO
425
426          DO  i = nxl, nxr
427             DO  k = nzb+1, nzt-1
428                w(k,nysg:nys-1,i) = nest_offl%w_south(0,k,i) + ddt_lsf * t_ref  *&
429                        ( nest_offl%w_south(1,k,i) - nest_offl%w_south(0,k,i) ) *&
430                           MERGE( 1.0_wp, 0.0_wp,                                &
431                                  BTEST( wall_flags_0(k,nysg:nys-1,i), 3 ) )
432             ENDDO
433          ENDDO
434
435          DO  i = nxlu, nxr
436             DO  k = nzb+1, nzt
437                u(k,nysg:nys-1,i) = nest_offl%u_south(0,k,i) + ddt_lsf * t_ref  *&
438                        ( nest_offl%u_south(1,k,i) - nest_offl%u_south(0,k,i) ) *&
439                           MERGE( 1.0_wp, 0.0_wp,                                &
440                                  BTEST( wall_flags_0(k,nysg:nys-1,i), 1 ) )
441             ENDDO
442          ENDDO
443
444          IF ( .NOT. neutral )  THEN
445             DO  i = nxl, nxr
446                DO  k = nzb+1, nzt
447                   pt(k,nysg:nys-1,i) = nest_offl%pt_south(0,k,i) + ddt_lsf *  &
448                                                                    t_ref   *  &
449                     ( nest_offl%pt_south(1,k,i) - nest_offl%pt_south(0,k,i) )
450 
451                ENDDO
452             ENDDO
453          ENDIF
454
455          IF ( humidity )  THEN
456             DO  i = nxl, nxr
457                DO  k = nzb+1, nzt
458                   q(k,nysg:nys-1,i) = nest_offl%q_south(0,k,i) + ddt_lsf   *  &
459                                                                    t_ref   *  &
460                       ( nest_offl%q_south(1,k,i) - nest_offl%q_south(0,k,i) )
461 
462                ENDDO
463             ENDDO
464          ENDIF
465
466       ENDIF
467
468       IF ( bc_dirichlet_n )  THEN
469
470          DO  i = nxl, nxr
471             DO  k = nzb+1, nzt
472                v(k,nyn+1:nyng,i)   = nest_offl%v_north(0,k,i) + ddt_lsf * t_ref *&
473                        ( nest_offl%v_north(1,k,i) - nest_offl%v_north(0,k,i) )  *&
474                           MERGE( 1.0_wp, 0.0_wp,                                 &
475                                  BTEST( wall_flags_0(k,nyn+1:nyng,i), 2 ) )
476             ENDDO
477          ENDDO
478          DO  i = nxl, nxr
479             DO  k = nzb+1, nzt-1
480                w(k,nyn+1:nyng,i) = nest_offl%w_north(0,k,i) + ddt_lsf * t_ref  *&
481                        ( nest_offl%w_north(1,k,i) - nest_offl%w_north(0,k,i) ) *&
482                           MERGE( 1.0_wp, 0.0_wp,                                &
483                                  BTEST( wall_flags_0(k,nyn+1:nyng,i), 3 ) )
484             ENDDO
485          ENDDO
486
487          DO  i = nxlu, nxr
488             DO  k = nzb+1, nzt
489                u(k,nyn+1:nyng,i) = nest_offl%u_north(0,k,i) + ddt_lsf * t_ref  *&
490                        ( nest_offl%u_north(1,k,i) - nest_offl%u_north(0,k,i) ) *&
491                           MERGE( 1.0_wp, 0.0_wp,                                &
492                                  BTEST( wall_flags_0(k,nyn+1:nyng,i), 1 ) )
493
494             ENDDO
495          ENDDO
496
497          IF ( .NOT. neutral )  THEN
498             DO  i = nxl, nxr
499                DO  k = nzb+1, nzt
500                   pt(k,nyn+1:nyng,i) = nest_offl%pt_north(0,k,i) + ddt_lsf *  &
501                                                                    t_ref   *  &
502                     ( nest_offl%pt_north(1,k,i) - nest_offl%pt_north(0,k,i) )
503 
504                ENDDO
505             ENDDO
506          ENDIF
507
508          IF ( humidity )  THEN
509             DO  i = nxl, nxr
510                DO  k = nzb+1, nzt
511                   q(k,nyn+1:nyng,i) = nest_offl%q_north(0,k,i) + ddt_lsf   *  &
512                                                                    t_ref   *  &
513                       ( nest_offl%q_north(1,k,i) - nest_offl%q_north(0,k,i) )
514 
515                ENDDO
516             ENDDO
517          ENDIF
518
519       ENDIF
520!
521!--    Top boundary.
522       DO  i = nxlu, nxr
523          DO  j = nys, nyn
524             u(nzt+1,j,i) = nest_offl%u_top(0,j,i) + ddt_lsf * t_ref *         &
525                        ( nest_offl%u_top(1,j,i) - nest_offl%u_top(0,j,i) ) *  &
526                           MERGE( 1.0_wp, 0.0_wp,                              &
527                                  BTEST( wall_flags_0(nzt+1,j,i), 1 ) )
528          ENDDO
529       ENDDO
530
531       DO  i = nxl, nxr
532          DO  j = nysv, nyn
533             v(nzt+1,j,i) = nest_offl%v_top(0,j,i) + ddt_lsf * t_ref *         &
534                        ( nest_offl%v_top(1,j,i) - nest_offl%v_top(0,j,i) ) *  &
535                           MERGE( 1.0_wp, 0.0_wp,                              &
536                                  BTEST( wall_flags_0(nzt+1,j,i), 2 ) )
537          ENDDO
538       ENDDO
539
540       DO  i = nxl, nxr
541          DO  j = nys, nyn
542             w(nzt:nzt+1,j,i) = nest_offl%w_top(0,j,i) + ddt_lsf * t_ref *     &
543                        ( nest_offl%w_top(1,j,i) - nest_offl%w_top(0,j,i) ) *  &
544                           MERGE( 1.0_wp, 0.0_wp,                              &
545                                  BTEST( wall_flags_0(nzt:nzt+1,j,i), 3 ) )
546          ENDDO
547       ENDDO
548
549
550       IF ( .NOT. neutral )  THEN
551          DO  i = nxl, nxr
552             DO  j = nys, nyn
553                pt(nzt+1,j,i) = nest_offl%pt_top(0,j,i) + ddt_lsf * t_ref *    &
554                        ( nest_offl%pt_top(1,j,i) - nest_offl%pt_top(0,j,i) )
555             ENDDO
556          ENDDO
557       ENDIF
558
559       IF ( humidity )  THEN
560          DO  i = nxl, nxr
561             DO  j = nys, nyn
562                q(nzt+1,j,i) = nest_offl%q_top(0,j,i) + ddt_lsf * t_ref *      &
563                        ( nest_offl%q_top(1,j,i) - nest_offl%q_top(0,j,i) )
564             ENDDO
565          ENDDO
566       ENDIF
567!
568!--    At the edges( left-south, left-north, right-south and right-north) set
569!--    data on ghost points.
570       IF ( bc_dirichlet_l  .AND.  bc_dirichlet_s )  THEN
571          DO  i = 1, nbgp
572             u(:,nys-i,nxlg:nxl)   = u(:,nys,nxlg:nxl)
573             w(:,nys-i,nxlg:nxl-1) = w(:,nys,nxlg:nxl-1)
574             IF ( .NOT. neutral )  pt(:,nys-i,nxlg:nxl-1) = pt(:,nys,nxlg:nxl-1)
575             IF ( humidity      )  q(:,nys-i,nxlg:nxl-1)  = q(:,nys,nxlg:nxl-1)
576          ENDDO
577          DO  i = 1, nbgp+1
578             v(:,nysv-i,nxlg:nxl-1) = v(:,nysv,nxlg:nxl-1)
579          ENDDO
580       ENDIF
581       IF ( bc_dirichlet_l  .AND.  bc_dirichlet_n )  THEN
582          DO  i = 1, nbgp
583             u(:,nyn+i,nxlg:nxl)   = u(:,nyn,nxlg:nxl)
584             v(:,nyn+i,nxlg:nxl-1) = v(:,nyn,nxlg:nxl-1)
585             w(:,nyn+i,nxlg:nxl-1) = w(:,nyn,nxlg:nxl-1)
586             IF ( .NOT. neutral )  pt(:,nyn+i,nxlg:nxl-1) = pt(:,nyn,nxlg:nxl-1)
587             IF ( humidity      )  q(:,nyn+i,nxlg:nxl-1)  = q(:,nyn,nxlg:nxl-1)
588          ENDDO
589       ENDIF
590       IF ( bc_dirichlet_r  .AND.  bc_dirichlet_s )  THEN
591          DO  i = 1, nbgp
592             u(:,nys-i,nxr+1:nxrg) = u(:,nys,nxr+1:nxrg)
593             w(:,nys-i,nxr+1:nxrg) = w(:,nys,nxr+1:nxrg)
594             IF ( .NOT. neutral )  pt(:,nys-i,nxr+1:nxrg) = pt(:,nys,nxr+1:nxrg)
595             IF ( humidity      )  q(:,nys-i,nxr+1:nxrg)  = q(:,nys,nxr+1:nxrg)
596          ENDDO
597          DO  i = 1, nbgp+1
598             v(:,nysv-i,nxr+1:nxrg) = v(:,nysv,nxr+1:nxrg)
599          ENDDO
600       ENDIF
601       IF ( bc_dirichlet_r  .AND.  bc_dirichlet_n )  THEN
602          DO  i = 1, nbgp
603             u(:,nyn+i,nxr+1:nxrg) = u(:,nyn,nxr+1:nxrg)
604             v(:,nyn+i,nxr+1:nxrg) = v(:,nyn,nxr+1:nxrg)
605             w(:,nyn+i,nxr+1:nxrg) = w(:,nyn,nxr+1:nxrg)
606             IF ( .NOT. neutral )  pt(:,nyn+i,nxr+1:nxrg) = pt(:,nyn,nxr+1:nxrg)
607             IF ( humidity      )  q(:,nyn+i,nxr+1:nxrg)  = q(:,nyn,nxr+1:nxrg)
608          ENDDO
609       ENDIF
610!
611!--    Moreover, set Neumann boundary condition for subgrid-scale TKE,
612!--    passive scalar, dissipation, and chemical species if required
613       IF ( rans_mode  .AND.  rans_tke_e )  THEN
614          IF (  bc_dirichlet_l )  diss(:,:,nxl-1) = diss(:,:,nxl)
615          IF (  bc_dirichlet_r )  diss(:,:,nxr+1) = diss(:,:,nxr)
616          IF (  bc_dirichlet_s )  diss(:,nys-1,:) = diss(:,nys,:)
617          IF (  bc_dirichlet_n )  diss(:,nyn+1,:) = diss(:,nyn,:)
618       ENDIF
619       IF ( .NOT. constant_diffusion )  THEN
620          IF (  bc_dirichlet_l )  e(:,:,nxl-1) = e(:,:,nxl)
621          IF (  bc_dirichlet_r )  e(:,:,nxr+1) = e(:,:,nxr)
622          IF (  bc_dirichlet_s )  e(:,nys-1,:) = e(:,nys,:)
623          IF (  bc_dirichlet_n )  e(:,nyn+1,:) = e(:,nyn,:)
624          e(nzt+1,:,:) = e(nzt,:,:)
625       ENDIF
626       IF ( passive_scalar )  THEN
627          IF (  bc_dirichlet_l )  s(:,:,nxl-1) = s(:,:,nxl)
628          IF (  bc_dirichlet_r )  s(:,:,nxr+1) = s(:,:,nxr)
629          IF (  bc_dirichlet_s )  s(:,nys-1,:) = s(:,nys,:)
630          IF (  bc_dirichlet_n )  s(:,nyn+1,:) = s(:,nyn,:)
631       ENDIF
632
633
634
635       CALL exchange_horiz( u, nbgp )
636       CALL exchange_horiz( v, nbgp )
637       CALL exchange_horiz( w, nbgp )
638       IF ( .NOT. neutral )  CALL exchange_horiz( pt, nbgp )
639       IF ( humidity      )  CALL exchange_horiz( q,  nbgp )
640
641!
642!--    Set surface pressure. Please note, time-dependent surface
643!--    pressure would require changes in anelastic approximation and
644!--    treatment of fluxes.
645!--    For the moment, comment this out!
646!      surface_pressure = nest_offl%surface_pressure(nest_offl%tind) +                 &
647!                                                      ddt_lsf * t_ref *       &
648!                                    ( nest_offl%surface_pressure(nest_offl%tind_p)    &
649!                                    - nest_offl%surface_pressure(nest_offl%tind) )
650
651    END SUBROUTINE lsf_nesting_offline
652
653!------------------------------------------------------------------------------!
654! Description:
655! ------------
656!> @todo Missing subroutine description.
657!------------------------------------------------------------------------------!
658    SUBROUTINE lsf_nudging_check_parameters
659
660       IMPLICIT NONE
661!
662!--    Check nudging and large scale forcing from external file
663       IF ( nudging  .AND.  (  .NOT.  large_scale_forcing ) )  THEN
664          message_string = 'Nudging requires large_scale_forcing = .T.. &'//   &
665                        'Surface fluxes and geostrophic wind should be &'//    &
666                        'prescribed in file LSF_DATA'
667          CALL message( 'check_parameters', 'PA0374', 1, 2, 0, 6, 0 )
668       ENDIF
669
670       IF ( large_scale_forcing  .AND.  ( bc_lr /= 'cyclic'  .OR.              &
671                                          bc_ns /= 'cyclic' ) )  THEN
672          message_string = 'Non-cyclic lateral boundaries do not allow for &'//&
673                        'the usage of large scale forcing from external file.'
674          CALL message( 'check_parameters', 'PA0375', 1, 2, 0, 6, 0 )
675       ENDIF
676
677       IF ( large_scale_forcing  .AND.  (  .NOT.  humidity ) )  THEN
678          message_string = 'The usage of large scale forcing from external &'//&
679                        'file LSF_DATA requires humidity = .T..'
680          CALL message( 'check_parameters', 'PA0376', 1, 2, 0, 6, 0 )
681       ENDIF
682
683       IF ( large_scale_forcing  .AND.  passive_scalar )  THEN
684          message_string = 'The usage of large scale forcing from external &'// &
685                        'file LSF_DATA is not implemented for passive scalars'
686          CALL message( 'check_parameters', 'PA0440', 1, 2, 0, 6, 0 )
687       ENDIF
688
689       IF ( large_scale_forcing  .AND.  topography /= 'flat'                   &
690                              .AND.  .NOT.  lsf_exception )  THEN
691          message_string = 'The usage of large scale forcing from external &'//&
692                        'file LSF_DATA is not implemented for non-flat topography'
693          CALL message( 'check_parameters', 'PA0377', 1, 2, 0, 6, 0 )
694       ENDIF
695
696       IF ( large_scale_forcing  .AND.  ocean )  THEN
697          message_string = 'The usage of large scale forcing from external &'//&
698                        'file LSF_DATA is not implemented for ocean runs'
699          CALL message( 'check_parameters', 'PA0378', 1, 2, 0, 6, 0 )
700       ENDIF
701
702    END SUBROUTINE lsf_nudging_check_parameters
703
704!------------------------------------------------------------------------------!
705! Description:
706! ------------
707!> Check data output of profiles for land surface model
708!------------------------------------------------------------------------------!
709    SUBROUTINE lsf_nudging_check_data_output_pr( variable, var_count, unit,    &
710                                                 dopr_unit )
711 
712       USE profil_parameter
713
714       IMPLICIT NONE
715   
716       CHARACTER (LEN=*) ::  unit      !<
717       CHARACTER (LEN=*) ::  variable  !<
718       CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
719 
720       INTEGER(iwp) ::  user_pr_index !<   
721       INTEGER(iwp) ::  var_count     !<
722
723       SELECT CASE ( TRIM( variable ) )
724       
725
726          CASE ( 'td_lsa_lpt' )
727             IF (  .NOT.  large_scale_forcing )  THEN
728                message_string = 'data_output_pr = ' //                        &
729                                 TRIM( data_output_pr(var_count) ) //          &
730                                 ' is not implemented for ' //                 &
731                                 'large_scale_forcing = .FALSE.'
732                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393',    &
733                               1, 2, 0, 6, 0 )
734             ELSE
735                dopr_index(var_count) = 81
736                dopr_unit             = 'K/s'
737                unit                  = 'K/s'
738                hom(:,2,81,:) = SPREAD( zu, 2, statistic_regions+1 )
739             ENDIF
740
741          CASE ( 'td_lsa_q' )
742             IF (  .NOT.  large_scale_forcing )  THEN
743                message_string = 'data_output_pr = ' //                        &
744                                 TRIM( data_output_pr(var_count) ) //          &
745                                 ' is not implemented for ' //                 &
746                                 'large_scale_forcing = .FALSE.'
747                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393',    &
748                               1, 2, 0, 6, 0 )
749             ELSE
750                dopr_index(var_count) = 82
751                dopr_unit             = 'kg/kgs'
752                unit                  = 'kg/kgs'
753                hom(:,2,82,:) = SPREAD( zu, 2, statistic_regions+1 )
754             ENDIF
755          CASE ( 'td_sub_lpt' )
756             IF (  .NOT.  large_scale_forcing )  THEN
757                message_string = 'data_output_pr = ' //                        &
758                                 TRIM( data_output_pr(var_count) ) //          &
759                                 ' is not implemented for ' //                 &
760                                 'large_scale_forcing = .FALSE.'
761                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393',    &
762                               1, 2, 0, 6, 0 )
763             ELSE
764                dopr_index(var_count) = 83
765                dopr_unit             = 'K/s'
766                unit                  = 'K/s'
767                hom(:,2,83,:) = SPREAD( zu, 2, statistic_regions+1 )
768             ENDIF
769
770          CASE ( 'td_sub_q' )
771             IF (  .NOT.  large_scale_forcing )  THEN
772                message_string = 'data_output_pr = ' //                        &
773                                 TRIM( data_output_pr(var_count) ) //          &
774                                 ' is not implemented for ' //                 &
775                                 'large_scale_forcing = .FALSE.'
776                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393',    &
777                               1, 2, 0, 6, 0 )
778             ELSE
779                dopr_index(var_count) = 84
780                dopr_unit             = 'kg/kgs'
781                unit                  = 'kg/kgs'
782                hom(:,2,84,:) = SPREAD( zu, 2, statistic_regions+1 )
783             ENDIF
784
785          CASE ( 'td_nud_lpt' )
786             IF (  .NOT.  nudging )  THEN
787                message_string = 'data_output_pr = ' //                        &
788                                 TRIM( data_output_pr(var_count) ) //          &
789                                 ' is not implemented for ' //                 &
790                                 'nudging = .FALSE.'
791                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0394',    &
792                               1, 2, 0, 6, 0 )
793             ELSE
794                dopr_index(var_count) = 85
795                dopr_unit             = 'K/s'
796                unit                  = 'K/s'
797                hom(:,2,85,:) = SPREAD( zu, 2, statistic_regions+1 )
798             ENDIF
799
800          CASE ( 'td_nud_q' )
801             IF (  .NOT.  nudging )  THEN
802                message_string = 'data_output_pr = ' //                        &
803                                 TRIM( data_output_pr(var_count) ) //          &
804                                 ' is not implemented for ' //                 &
805                                 'nudging = .FALSE.'
806                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0394',    &
807                               1, 2, 0, 6, 0 )
808             ELSE
809                dopr_index(var_count) = 86
810                dopr_unit             = 'kg/kgs'
811                unit                  = 'kg/kgs'
812                hom(:,2,86,:) = SPREAD( zu, 2, statistic_regions+1 )
813             ENDIF
814
815          CASE ( 'td_nud_u' )
816             IF (  .NOT.  nudging )  THEN
817                message_string = 'data_output_pr = ' //                        &
818                                 TRIM( data_output_pr(var_count) ) //          &
819                                 ' is not implemented for ' //                 &
820                                 'nudging = .FALSE.'
821                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0394',    &
822                               1, 2, 0, 6, 0 )
823             ELSE
824                dopr_index(var_count) = 87
825                dopr_unit             = 'm/s2'
826                unit                  = 'm/s2'
827                hom(:,2,87,:) = SPREAD( zu, 2, statistic_regions+1 )
828             ENDIF
829
830          CASE ( 'td_nud_v' )
831             IF (  .NOT.  nudging )  THEN
832                message_string = 'data_output_pr = ' //                        &
833                                 TRIM( data_output_pr(var_count) ) //          &
834                                 ' is not implemented for ' //                 &
835                                 'nudging = .FALSE.'
836                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0394',    &
837                               1, 2, 0, 6, 0 )
838             ELSE
839                dopr_index(var_count) = 88
840                dopr_unit             = 'm/s2'
841                unit                  = 'm/s2'
842                hom(:,2,88,:) = SPREAD( zu, 2, statistic_regions+1 )
843             ENDIF
844
845
846          CASE DEFAULT
847             unit = 'illegal'
848   
849       END SELECT
850
851    END SUBROUTINE lsf_nudging_check_data_output_pr
852
853!------------------------------------------------------------------------------!
854! Description:
855! ------------
856!> @todo Missing subroutine description.
857!------------------------------------------------------------------------------!
858    SUBROUTINE lsf_nudging_header ( io )
859
860       IMPLICIT NONE
861
862       INTEGER(iwp), INTENT(IN) ::  io !< Unit of the output file
863
864       WRITE ( io, 1 )
865       IF ( large_scale_forcing )  THEN
866          WRITE ( io, 3 )
867          WRITE ( io, 4 )
868
869          IF ( large_scale_subsidence )  THEN
870             IF ( .NOT. use_subsidence_tendencies )  THEN
871                WRITE ( io, 5 )
872             ELSE
873                WRITE ( io, 6 )
874             ENDIF
875          ENDIF
876
877          IF ( bc_pt_b == 'dirichlet' )  THEN
878             WRITE ( io, 12 )
879          ELSEIF ( bc_pt_b == 'neumann' )  THEN
880             WRITE ( io, 13 )
881          ENDIF
882
883          IF ( bc_q_b == 'dirichlet' )  THEN
884             WRITE ( io, 14 )
885          ELSEIF ( bc_q_b == 'neumann' )  THEN
886             WRITE ( io, 15 )
887          ENDIF
888
889          WRITE ( io, 7 )
890          IF ( nudging )  THEN
891             WRITE ( io, 10 )
892          ENDIF
893       ELSE
894          WRITE ( io, 2 )
895          WRITE ( io, 11 )
896       ENDIF
897       IF ( large_scale_subsidence )  THEN
898          WRITE ( io, 8 )
899          WRITE ( io, 9 )
900       ENDIF
901
902
9031 FORMAT (//' Large scale forcing and nudging:'/ &
904              ' -------------------------------'/)
9052 FORMAT (' --> No large scale forcing from external is used (default) ')
9063 FORMAT (' --> Large scale forcing from external file LSF_DATA is used: ')
9074 FORMAT ('     - large scale advection tendencies ')
9085 FORMAT ('     - large scale subsidence velocity w_subs ')
9096 FORMAT ('     - large scale subsidence tendencies ')
9107 FORMAT ('     - and geostrophic wind components ug and vg')
9118 FORMAT (' --> Large-scale vertical motion is used in the ', &
912                  'prognostic equation(s) for')
9139 FORMAT ('     the scalar(s) only')
91410 FORMAT (' --> Nudging is used')
91511 FORMAT (' --> No nudging is used (default) ')
91612 FORMAT ('     - prescribed surface values for temperature')
91713 FORMAT ('     - prescribed surface fluxes for temperature')
91814 FORMAT ('     - prescribed surface values for humidity')
91915 FORMAT ('     - prescribed surface fluxes for humidity')
920
921    END SUBROUTINE lsf_nudging_header
922
923!------------------------------------------------------------------------------!
924! Description:
925! ------------
926!> @todo Missing subroutine description.
927!------------------------------------------------------------------------------!
928    SUBROUTINE lsf_init
929
930       USE control_parameters,                                                 &
931           ONLY:  bc_lr_cyc, bc_ns_cyc
932
933       USE netcdf_data_input_mod,                                              &
934           ONLY:  netcdf_data_input_lsf 
935
936       IMPLICIT NONE
937
938       CHARACTER(100) ::  chmess      !<
939       CHARACTER(1)   ::  hash        !<
940
941       INTEGER(iwp) ::  ierrn         !<
942       INTEGER(iwp) ::  finput = 90   !<
943       INTEGER(iwp) ::  k             !<
944       INTEGER(iwp) ::  nt            !<
945       INTEGER(iwp) ::  t             !< running index for time levels
946
947       REAL(wp) ::  fac               !<
948       REAL(wp) ::  highheight        !<
949       REAL(wp) ::  highug_vert       !<
950       REAL(wp) ::  highvg_vert       !<
951       REAL(wp) ::  highwsubs_vert    !<
952       REAL(wp) ::  lowheight         !<
953       REAL(wp) ::  lowug_vert        !<
954       REAL(wp) ::  lowvg_vert        !<
955       REAL(wp) ::  lowwsubs_vert     !<
956       REAL(wp) ::  high_td_lsa_lpt   !<
957       REAL(wp) ::  low_td_lsa_lpt    !<
958       REAL(wp) ::  high_td_lsa_q     !<
959       REAL(wp) ::  low_td_lsa_q      !<
960       REAL(wp) ::  high_td_sub_lpt   !<
961       REAL(wp) ::  low_td_sub_lpt    !<
962       REAL(wp) ::  high_td_sub_q     !<
963       REAL(wp) ::  low_td_sub_q      !<
964       REAL(wp) ::  r_dummy           !<
965
966       IF ( nesting_offline )  THEN
967!
968!--       Allocate arrays for geostrophic wind components. Arrays will
969!--       incorporate 2 time levels in order to interpolate in between. Please
970!--       note, forcing using geostrophic wind components is only required in
971!--       case of cyclic boundary conditions.
972          IF ( bc_lr_cyc  .AND.  bc_ns_cyc )  THEN
973             ALLOCATE( nest_offl%ug(0:1,nzb:nzt+1) )
974             ALLOCATE( nest_offl%vg(0:1,nzb:nzt+1) )
975          ENDIF
976!
977!--       Allocate arrays for reading boundary values. Arrays will incorporate 2
978!--       time levels in order to interpolate in between.
979          IF ( bc_dirichlet_l )  THEN
980             ALLOCATE( nest_offl%u_left(0:1,nzb+1:nzt,nys:nyn)  )
981             ALLOCATE( nest_offl%v_left(0:1,nzb+1:nzt,nysv:nyn) )
982             ALLOCATE( nest_offl%w_left(0:1,nzb+1:nzt-1,nys:nyn) )
983             IF ( humidity )      ALLOCATE( nest_offl%q_left(0:1,nzb+1:nzt,nys:nyn)  )
984             IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_left(0:1,nzb+1:nzt,nys:nyn) )
985          ENDIF
986          IF ( bc_dirichlet_r )  THEN
987             ALLOCATE( nest_offl%u_right(0:1,nzb+1:nzt,nys:nyn)  )
988             ALLOCATE( nest_offl%v_right(0:1,nzb+1:nzt,nysv:nyn) )
989             ALLOCATE( nest_offl%w_right(0:1,nzb+1:nzt-1,nys:nyn) )
990             IF ( humidity )      ALLOCATE( nest_offl%q_right(0:1,nzb+1:nzt,nys:nyn)  )
991             IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_right(0:1,nzb+1:nzt,nys:nyn) )
992          ENDIF
993          IF ( bc_dirichlet_n )  THEN
994             ALLOCATE( nest_offl%u_north(0:1,nzb+1:nzt,nxlu:nxr) )
995             ALLOCATE( nest_offl%v_north(0:1,nzb+1:nzt,nxl:nxr)  )
996             ALLOCATE( nest_offl%w_north(0:1,nzb+1:nzt-1,nxl:nxr) )
997             IF ( humidity )      ALLOCATE( nest_offl%q_north(0:1,nzb+1:nzt,nxl:nxr)  )
998             IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_north(0:1,nzb+1:nzt,nxl:nxr) )
999          ENDIF
1000          IF ( bc_dirichlet_s )  THEN
1001             ALLOCATE( nest_offl%u_south(0:1,nzb+1:nzt,nxlu:nxr) )
1002             ALLOCATE( nest_offl%v_south(0:1,nzb+1:nzt,nxl:nxr)  )
1003             ALLOCATE( nest_offl%w_south(0:1,nzb+1:nzt-1,nxl:nxr)    )
1004             IF ( humidity )      ALLOCATE( nest_offl%q_south(0:1,nzb+1:nzt,nxl:nxr)  )
1005             IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_south(0:1,nzb+1:nzt,nxl:nxr) )
1006          ENDIF
1007         
1008          ALLOCATE( nest_offl%u_top(0:1,nys:nyn,nxlu:nxr) )
1009          ALLOCATE( nest_offl%v_top(0:1,nysv:nyn,nxl:nxr) )
1010          ALLOCATE( nest_offl%w_top(0:1,nys:nyn,nxl:nxr)  )
1011          IF ( humidity )      ALLOCATE( nest_offl%q_top(0:1,nys:nyn,nxl:nxr)  )
1012          IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_top(0:1,nys:nyn,nxl:nxr) )
1013
1014!
1015!--       Read COSMO data at lateral and top boundaries
1016          CALL netcdf_data_input_lsf
1017!
1018!--       Write COSMO data at lateral and top boundaries
1019          CALL lsf_nesting_offline
1020!
1021!--       After 3D data is initialized, ensure mass conservation
1022          CALL lsf_nesting_offline_mass_conservation
1023!
1024!--       Initialize surface pressure. Please note, time-dependent surface
1025!--       pressure would require changes in anelastic approximation and
1026!--       treatment of fluxes.
1027!--       For the moment, comment this out!
1028!         surface_pressure = nest_offl%surface_pressure(0)
1029
1030       ELSE
1031
1032          ALLOCATE( p_surf(0:nlsf), pt_surf(0:nlsf), q_surf(0:nlsf),           &
1033                    qsws_surf(0:nlsf), shf_surf(0:nlsf),                       &
1034                    td_lsa_lpt(nzb:nzt+1,0:nlsf), td_lsa_q(nzb:nzt+1,0:nlsf),  &
1035                    td_sub_lpt(nzb:nzt+1,0:nlsf), td_sub_q(nzb:nzt+1,0:nlsf),  &
1036                    time_vert(0:nlsf), time_surf(0:nlsf),                      &
1037                    ug_vert(nzb:nzt+1,0:nlsf), vg_vert(nzb:nzt+1,0:nlsf),      &
1038                    wsubs_vert(nzb:nzt+1,0:nlsf) )
1039
1040          p_surf = 0.0_wp; pt_surf = 0.0_wp; q_surf = 0.0_wp; qsws_surf = 0.0_wp
1041          shf_surf = 0.0_wp; time_vert = 0.0_wp; td_lsa_lpt = 0.0_wp
1042          td_lsa_q = 0.0_wp; td_sub_lpt = 0.0_wp; td_sub_q = 0.0_wp
1043          time_surf = 0.0_wp; ug_vert = 0.0_wp; vg_vert = 0.0_wp
1044          wsubs_vert = 0.0_wp
1045
1046!
1047!--       Array for storing large scale forcing and nudging tendencies at each
1048!--       timestep for data output
1049          ALLOCATE( sums_ls_l(nzb:nzt+1,0:7) )
1050          sums_ls_l = 0.0_wp
1051
1052          ngp_sums_ls = (nz+2)*6
1053
1054          OPEN ( finput, FILE='LSF_DATA', STATUS='OLD', &
1055                 FORM='FORMATTED', IOSTAT=ierrn )
1056
1057          IF ( ierrn /= 0 )  THEN
1058             message_string = 'file LSF_DATA does not exist'
1059             CALL message( 'ls_forcing', 'PA0368', 1, 2, 0, 6, 0 )
1060          ENDIF
1061
1062          ierrn = 0
1063!
1064!--       First three lines of LSF_DATA contain header
1065          READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
1066          READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
1067          READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
1068
1069          IF ( ierrn /= 0 )  THEN
1070             message_string = 'errors in file LSF_DATA'
1071             CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
1072          ENDIF
1073
1074!
1075!--       Surface values are read in
1076          nt     = 0
1077          ierrn = 0
1078
1079          DO WHILE ( time_surf(nt) < end_time )
1080             nt = nt + 1
1081             READ ( finput, *, IOSTAT = ierrn ) time_surf(nt), shf_surf(nt),   &
1082                                                qsws_surf(nt), pt_surf(nt),    &
1083                                                q_surf(nt), p_surf(nt)
1084
1085             IF ( ierrn /= 0 )  THEN
1086               WRITE ( message_string, * ) 'No time dependent surface ' //     &
1087                                 'variables in & LSF_DATA for end of run found'
1088
1089                CALL message( 'ls_forcing', 'PA0363', 1, 2, 0, 6, 0 )
1090             ENDIF
1091          ENDDO
1092
1093          IF ( time_surf(1) > end_time )  THEN
1094             WRITE ( message_string, * ) 'Time dependent surface variables in ' // &
1095                                         '&LSF_DATA set in after end of ' ,        &
1096                                         'simulation - lsf_surf is set to FALSE'
1097             CALL message( 'ls_forcing', 'PA0371', 0, 0, 0, 6, 0 )
1098             lsf_surf = .FALSE.
1099          ENDIF
1100
1101!
1102!--       Go to the end of the list with surface variables
1103          DO WHILE ( ierrn == 0 )
1104             READ ( finput, *, IOSTAT = ierrn ) r_dummy
1105          ENDDO
1106
1107!
1108!--       Profiles of ug, vg and w_subs are read in (large scale forcing)
1109
1110          nt = 0
1111          DO WHILE ( time_vert(nt) < end_time )
1112             nt = nt + 1
1113             hash = "#"
1114             ierrn = 1 ! not zero
1115!
1116!--          Search for the next line consisting of "# time",
1117!--          from there onwards the profiles will be read
1118             DO WHILE ( .NOT. ( hash == "#" .AND. ierrn == 0 ) ) 
1119                READ ( finput, *, IOSTAT=ierrn ) hash, time_vert(nt)
1120                IF ( ierrn < 0 )  THEN
1121                   WRITE( message_string, * ) 'No time dependent vertical profiles',&
1122                                    ' in & LSF_DATA for end of run found'
1123                   CALL message( 'ls_forcing', 'PA0372', 1, 2, 0, 6, 0 )
1124                ENDIF
1125             ENDDO
1126
1127             IF ( nt == 1 .AND. time_vert(nt) > end_time ) EXIT
1128
1129             READ ( finput, *, IOSTAT=ierrn ) lowheight, lowug_vert, lowvg_vert,&
1130                                              lowwsubs_vert, low_td_lsa_lpt,    &
1131                                              low_td_lsa_q, low_td_sub_lpt,     &
1132                                              low_td_sub_q
1133             IF ( ierrn /= 0 )  THEN
1134                message_string = 'errors in file LSF_DATA'
1135                CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
1136             ENDIF
1137
1138             READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert,         &
1139                                              highvg_vert, highwsubs_vert,     &
1140                                              high_td_lsa_lpt, high_td_lsa_q,  &
1141                                              high_td_sub_lpt, high_td_sub_q
1142         
1143             IF ( ierrn /= 0 )  THEN
1144                message_string = 'errors in file LSF_DATA'
1145                CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
1146             ENDIF
1147
1148
1149             DO  k = nzb, nzt+1
1150                IF ( highheight < zu(k) )  THEN
1151                   lowheight      = highheight
1152                   lowug_vert     = highug_vert
1153                   lowvg_vert     = highvg_vert
1154                   lowwsubs_vert  = highwsubs_vert
1155                   low_td_lsa_lpt = high_td_lsa_lpt
1156                   low_td_lsa_q   = high_td_lsa_q
1157                   low_td_sub_lpt = high_td_sub_lpt
1158                   low_td_sub_q   = high_td_sub_q
1159
1160                   ierrn = 0
1161                   READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert,    &
1162                                                    highvg_vert, highwsubs_vert,&
1163                                                    high_td_lsa_lpt,            &
1164                                                    high_td_lsa_q,              &
1165                                                    high_td_sub_lpt, high_td_sub_q
1166
1167                   IF ( ierrn /= 0 )  THEN
1168                      WRITE( message_string, * ) 'zu(',k,') = ', zu(k), 'm ',  &
1169                           'is higher than the maximum height in LSF_DATA ',   &
1170                           'which is ', lowheight, 'm. Interpolation on PALM ',&
1171                           'grid is not possible.'
1172                      CALL message( 'ls_forcing', 'PA0395', 1, 2, 0, 6, 0 )
1173                   ENDIF
1174
1175                ENDIF
1176
1177!
1178!--             Interpolation of prescribed profiles in space
1179                fac = (highheight-zu(k))/(highheight - lowheight)
1180
1181                ug_vert(k,nt)    = fac * lowug_vert                            &
1182                                   + ( 1.0_wp - fac ) * highug_vert
1183                vg_vert(k,nt)    = fac * lowvg_vert                            &
1184                                   + ( 1.0_wp - fac ) * highvg_vert
1185                wsubs_vert(k,nt) = fac * lowwsubs_vert                         &
1186                                   + ( 1.0_wp - fac ) * highwsubs_vert
1187
1188                td_lsa_lpt(k,nt) = fac * low_td_lsa_lpt                        &
1189                                   + ( 1.0_wp - fac ) * high_td_lsa_lpt
1190                td_lsa_q(k,nt)   = fac * low_td_lsa_q                          &
1191                                   + ( 1.0_wp - fac ) * high_td_lsa_q
1192                td_sub_lpt(k,nt) = fac * low_td_sub_lpt                        &
1193                                   + ( 1.0_wp - fac ) * high_td_sub_lpt
1194                td_sub_q(k,nt)   = fac * low_td_sub_q                          &
1195                                   + ( 1.0_wp - fac ) * high_td_sub_q
1196
1197             ENDDO
1198
1199          ENDDO
1200
1201!
1202!--       Large scale vertical velocity has to be zero at the surface
1203          wsubs_vert(nzb,:) = 0.0_wp
1204   
1205          IF ( time_vert(1) > end_time )  THEN
1206             WRITE ( message_string, * ) 'Time dependent large scale profile ',&
1207                                'forcing from&LSF_DATA sets in after end of ' ,&
1208                                'simulation - lsf_vert is set to FALSE'
1209             CALL message( 'ls_forcing', 'PA0373', 0, 0, 0, 6, 0 )
1210             lsf_vert = .FALSE.
1211          ENDIF
1212
1213          CLOSE( finput )
1214
1215       ENDIF
1216
1217    END SUBROUTINE lsf_init
1218
1219!------------------------------------------------------------------------------!
1220! Description:
1221! ------------
1222!> @todo Missing subroutine description.
1223!------------------------------------------------------------------------------!
1224    SUBROUTINE ls_forcing_surf ( time )
1225
1226       IMPLICIT NONE
1227
1228       INTEGER(iwp) ::  nt                     !<
1229
1230       REAL(wp)             :: dum_surf_flux  !<
1231       REAL(wp)             :: fac            !<
1232       REAL(wp), INTENT(in) :: time           !<
1233
1234!
1235!--    Interpolation in time of LSF_DATA at the surface
1236       nt = 1
1237       DO WHILE ( time > time_surf(nt) )
1238          nt = nt + 1
1239       ENDDO
1240       IF ( time /= time_surf(nt) )  THEN
1241         nt = nt - 1
1242       ENDIF
1243
1244       fac = ( time -time_surf(nt) ) / ( time_surf(nt+1) - time_surf(nt) )
1245
1246       IF ( ibc_pt_b == 0 )  THEN
1247!
1248!--       In case of Dirichlet boundary condition shf must not
1249!--       be set - it is calculated via MOST in prandtl_fluxes
1250          pt_surface = pt_surf(nt) + fac * ( pt_surf(nt+1) - pt_surf(nt) )
1251
1252       ELSEIF ( ibc_pt_b == 1 )  THEN
1253!
1254!--       In case of Neumann boundary condition pt_surface is needed for
1255!--       calculation of reference density
1256          dum_surf_flux = ( shf_surf(nt) + fac *                               &
1257                            ( shf_surf(nt+1) - shf_surf(nt) )                  &
1258                          ) * heatflux_input_conversion(nzb)
1259!
1260!--       Save surface sensible heat flux on default, natural and urban surface
1261!--       type, if required
1262          IF ( surf_def_h(0)%ns >= 1 )  surf_def_h(0)%shf(:) = dum_surf_flux
1263          IF ( surf_lsm_h%ns    >= 1 )  surf_lsm_h%shf(:)    = dum_surf_flux
1264          IF ( surf_usm_h%ns    >= 1 )  surf_usm_h%shf(:)    = dum_surf_flux
1265
1266          pt_surface    = pt_surf(nt) + fac * ( pt_surf(nt+1) - pt_surf(nt) )
1267
1268       ENDIF
1269
1270       IF ( ibc_q_b == 0 )  THEN
1271!
1272!--       In case of Dirichlet boundary condition qsws must not
1273!--       be set - it is calculated via MOST in prandtl_fluxes
1274          q_surface = q_surf(nt) + fac * ( q_surf(nt+1) - q_surf(nt) )
1275
1276       ELSEIF ( ibc_q_b == 1 )  THEN
1277          dum_surf_flux = ( qsws_surf(nt) + fac *                              &
1278                             ( qsws_surf(nt+1) - qsws_surf(nt) )               &
1279                             ) * waterflux_input_conversion(nzb)
1280!
1281!--       Save surface sensible heat flux on default, natural and urban surface
1282!--       type, if required
1283          IF ( surf_def_h(0)%ns >= 1 )  surf_def_h(0)%qsws(:) = dum_surf_flux
1284          IF ( surf_lsm_h%ns    >= 1 )  surf_lsm_h%qsws(:)    = dum_surf_flux
1285          IF ( surf_usm_h%ns    >= 1 )  surf_usm_h%qsws(:)    = dum_surf_flux
1286
1287       ENDIF
1288!
1289!--    Surface heat- and waterflux will be written later onto surface elements
1290       IF ( .NOT.  neutral  .AND.  constant_heatflux  .AND.                    &
1291            TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1292             surface_heatflux = shf_surf(1)
1293       ENDIF
1294       IF ( humidity  .AND.  constant_waterflux  .AND.                         &
1295            TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1296             surface_waterflux = qsws_surf(1)
1297       ENDIF
1298
1299       surface_pressure = p_surf(nt) + fac * ( p_surf(nt+1) - p_surf(nt) )
1300
1301    END SUBROUTINE ls_forcing_surf
1302
1303
1304
1305
1306!------------------------------------------------------------------------------!
1307! Description:
1308! ------------
1309!> @todo Missing subroutine description.
1310!------------------------------------------------------------------------------!
1311    SUBROUTINE ls_forcing_vert ( time )
1312
1313
1314       IMPLICIT NONE
1315
1316       INTEGER(iwp) ::  nt                     !<
1317
1318       REAL(wp)             ::  fac           !<
1319       REAL(wp), INTENT(in) ::  time          !<
1320
1321!
1322!--    Interpolation in time of LSF_DATA for ug, vg and w_subs
1323       nt = 1
1324       DO WHILE ( time > time_vert(nt) )
1325          nt = nt + 1
1326       ENDDO
1327       IF ( time /= time_vert(nt) )  THEN
1328         nt = nt - 1
1329       ENDIF
1330
1331       fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
1332
1333       ug     = ug_vert(:,nt) + fac * ( ug_vert(:,nt+1) - ug_vert(:,nt) )
1334       vg     = vg_vert(:,nt) + fac * ( vg_vert(:,nt+1) - vg_vert(:,nt) )
1335
1336       IF ( large_scale_subsidence )  THEN
1337          w_subs = wsubs_vert(:,nt)                                            &
1338                   + fac * ( wsubs_vert(:,nt+1) - wsubs_vert(:,nt) )
1339       ENDIF
1340
1341    END SUBROUTINE ls_forcing_vert
1342
1343
1344!------------------------------------------------------------------------------!
1345! Description:
1346! ------------
1347!> Call for all grid points
1348!------------------------------------------------------------------------------!
1349    SUBROUTINE ls_advec ( time, prog_var )
1350     
1351
1352       IMPLICIT NONE
1353
1354       CHARACTER (LEN=*) ::  prog_var   !<
1355
1356       REAL(wp), INTENT(in)  :: time    !<
1357       REAL(wp) :: fac                  !< 
1358
1359       INTEGER(iwp) ::  i               !<
1360       INTEGER(iwp) ::  j               !<
1361       INTEGER(iwp) ::  k               !<
1362       INTEGER(iwp) ::  nt               !<
1363
1364!
1365!--    Interpolation in time of LSF_DATA
1366       nt = 1
1367       DO WHILE ( time > time_vert(nt) )
1368          nt = nt + 1
1369       ENDDO
1370       IF ( time /= time_vert(nt) )  THEN
1371         nt = nt - 1
1372       ENDIF
1373
1374       fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
1375
1376!
1377!--    Add horizontal large scale advection tendencies of pt and q
1378       SELECT CASE ( prog_var )
1379
1380          CASE ( 'pt' )
1381
1382             DO  i = nxl, nxr
1383                DO  j = nys, nyn
1384                   DO  k = nzb+1, nzt
1385                      tend(k,j,i) = tend(k,j,i) + td_lsa_lpt(k,nt) + fac *     &
1386                                    ( td_lsa_lpt(k,nt+1) - td_lsa_lpt(k,nt) ) *&
1387                                        MERGE( 1.0_wp, 0.0_wp,                 &
1388                                               BTEST( wall_flags_0(k,j,i), 0 ) )
1389                   ENDDO
1390                ENDDO
1391             ENDDO
1392
1393          CASE ( 'q' )
1394
1395             DO  i = nxl, nxr
1396                DO  j = nys, nyn
1397                   DO  k = nzb+1, nzt
1398                      tend(k,j,i) = tend(k,j,i) + td_lsa_q(k,nt) + fac *       &
1399                                    ( td_lsa_q(k,nt+1) - td_lsa_q(k,nt) ) *    &
1400                                        MERGE( 1.0_wp, 0.0_wp,                 &
1401                                               BTEST( wall_flags_0(k,j,i), 0 ) )
1402                   ENDDO
1403                ENDDO
1404             ENDDO
1405
1406       END SELECT
1407
1408!
1409!--    Subsidence of pt and q with prescribed subsidence tendencies
1410       IF ( large_scale_subsidence .AND. use_subsidence_tendencies )  THEN
1411
1412          SELECT CASE ( prog_var )
1413
1414             CASE ( 'pt' )
1415
1416                DO  i = nxl, nxr
1417                   DO  j = nys, nyn
1418                      DO  k = nzb+1, nzt
1419                         tend(k,j,i) = tend(k,j,i) + td_sub_lpt(k,nt) + fac *  &
1420                                     ( td_sub_lpt(k,nt+1) - td_sub_lpt(k,nt) )*&
1421                                        MERGE( 1.0_wp, 0.0_wp,                 &
1422                                               BTEST( wall_flags_0(k,j,i), 0 ) )
1423                      ENDDO
1424                   ENDDO
1425                ENDDO
1426 
1427             CASE ( 'q' )
1428
1429                DO  i = nxl, nxr
1430                   DO  j = nys, nyn
1431                      DO  k = nzb+1, nzt
1432                         tend(k,j,i) = tend(k,j,i) + td_sub_q(k,nt) + fac *    &
1433                                       ( td_sub_q(k,nt+1) - td_sub_q(k,nt) ) * &
1434                                        MERGE( 1.0_wp, 0.0_wp,                 &
1435                                               BTEST( wall_flags_0(k,j,i), 0 ) )
1436                      ENDDO
1437                   ENDDO
1438                ENDDO
1439
1440          END SELECT
1441
1442       ENDIF
1443
1444    END SUBROUTINE ls_advec
1445
1446
1447!------------------------------------------------------------------------------!
1448! Description:
1449! ------------
1450!> Call for grid point i,j
1451!------------------------------------------------------------------------------!
1452    SUBROUTINE ls_advec_ij ( i, j, time, prog_var )
1453
1454       IMPLICIT NONE
1455
1456       CHARACTER (LEN=*) ::  prog_var   !<
1457
1458       REAL(wp), INTENT(in)  :: time    !<
1459       REAL(wp) :: fac                  !<
1460
1461       INTEGER(iwp) ::  i               !<
1462       INTEGER(iwp) ::  j               !<
1463       INTEGER(iwp) ::  k               !<
1464       INTEGER(iwp) ::  nt               !<
1465
1466!
1467!--    Interpolation in time of LSF_DATA
1468       nt = 1
1469       DO WHILE ( time > time_vert(nt) )
1470          nt = nt + 1
1471       ENDDO
1472       IF ( time /= time_vert(nt) )  THEN
1473         nt = nt - 1
1474       ENDIF
1475
1476       fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
1477
1478!
1479!--    Add horizontal large scale advection tendencies of pt and q
1480       SELECT CASE ( prog_var )
1481
1482          CASE ( 'pt' )
1483
1484             DO  k = nzb+1, nzt
1485                tend(k,j,i) = tend(k,j,i) + td_lsa_lpt(k,nt)                   &
1486                             + fac * ( td_lsa_lpt(k,nt+1) - td_lsa_lpt(k,nt) )*&
1487                                        MERGE( 1.0_wp, 0.0_wp,                 &
1488                                               BTEST( wall_flags_0(k,j,i), 0 ) )
1489             ENDDO
1490
1491          CASE ( 'q' )
1492
1493             DO  k = nzb+1, nzt
1494                tend(k,j,i) = tend(k,j,i) + td_lsa_q(k,nt)                     &
1495                              + fac * ( td_lsa_q(k,nt+1) - td_lsa_q(k,nt) ) *  &
1496                                        MERGE( 1.0_wp, 0.0_wp,                 &
1497                                               BTEST( wall_flags_0(k,j,i), 0 ) )
1498             ENDDO
1499
1500       END SELECT
1501
1502!
1503!--    Subsidence of pt and q with prescribed profiles
1504       IF ( large_scale_subsidence .AND. use_subsidence_tendencies )  THEN
1505
1506          SELECT CASE ( prog_var )
1507
1508             CASE ( 'pt' )
1509
1510                DO  k = nzb+1, nzt
1511                   tend(k,j,i) = tend(k,j,i) + td_sub_lpt(k,nt) + fac *        &
1512                                 ( td_sub_lpt(k,nt+1) - td_sub_lpt(k,nt) ) *   &
1513                                        MERGE( 1.0_wp, 0.0_wp,                 &
1514                                               BTEST( wall_flags_0(k,j,i), 0 ) )
1515                ENDDO
1516 
1517             CASE ( 'q' )
1518
1519                DO  k = nzb+1, nzt
1520                   tend(k,j,i) = tend(k,j,i) + td_sub_q(k,nt) + fac *          &
1521                                 ( td_sub_q(k,nt+1) - td_sub_q(k,nt) ) *       &
1522                                        MERGE( 1.0_wp, 0.0_wp,                 &
1523                                               BTEST( wall_flags_0(k,j,i), 0 ) )
1524                ENDDO
1525
1526          END SELECT
1527
1528       ENDIF
1529
1530    END SUBROUTINE ls_advec_ij
1531
1532
1533!------------------------------------------------------------------------------!
1534! Description:
1535! ------------
1536!> @todo Missing subroutine description.
1537!------------------------------------------------------------------------------!
1538    SUBROUTINE nudge_init
1539
1540       IMPLICIT NONE
1541
1542
1543       INTEGER(iwp) ::  finput = 90  !<
1544       INTEGER(iwp) ::  ierrn        !<
1545       INTEGER(iwp) ::  k            !<
1546       INTEGER(iwp) ::  nt            !<
1547
1548       CHARACTER(1) ::  hash     !<
1549
1550       REAL(wp) ::  highheight   !<
1551       REAL(wp) ::  highqnudge   !<
1552       REAL(wp) ::  highptnudge  !<
1553       REAL(wp) ::  highunudge   !<
1554       REAL(wp) ::  highvnudge   !<
1555       REAL(wp) ::  highwnudge   !<
1556       REAL(wp) ::  hightnudge   !<
1557
1558       REAL(wp) ::  lowheight    !<
1559       REAL(wp) ::  lowqnudge    !<
1560       REAL(wp) ::  lowptnudge   !<
1561       REAL(wp) ::  lowunudge    !<
1562       REAL(wp) ::  lowvnudge    !<
1563       REAL(wp) ::  lowwnudge    !<
1564       REAL(wp) ::  lowtnudge    !<
1565
1566       REAL(wp) ::  fac          !<
1567
1568       ALLOCATE( ptnudge(nzb:nzt+1,1:ntnudge), qnudge(nzb:nzt+1,1:ntnudge), &
1569                 tnudge(nzb:nzt+1,1:ntnudge), unudge(nzb:nzt+1,1:ntnudge),  &
1570                 vnudge(nzb:nzt+1,1:ntnudge), wnudge(nzb:nzt+1,1:ntnudge)  )
1571
1572       ALLOCATE( tmp_tnudge(nzb:nzt) )
1573
1574       ALLOCATE( timenudge(0:ntnudge) )
1575
1576       ptnudge = 0.0_wp; qnudge = 0.0_wp; tnudge = 0.0_wp; unudge = 0.0_wp
1577       vnudge = 0.0_wp; wnudge = 0.0_wp; timenudge = 0.0_wp
1578!
1579!--    Initialize array tmp_nudge with a current nudging time scale of 6 hours
1580       tmp_tnudge = 21600.0_wp
1581
1582       nt = 0
1583       OPEN ( finput, FILE='NUDGING_DATA', STATUS='OLD', &
1584              FORM='FORMATTED', IOSTAT=ierrn )
1585
1586       IF ( ierrn /= 0 )  THEN
1587          message_string = 'file NUDGING_DATA does not exist'
1588          CALL message( 'nudging', 'PA0365', 1, 2, 0, 6, 0 )
1589       ENDIF
1590
1591       ierrn = 0
1592
1593 rloop:DO
1594          nt = nt + 1
1595          hash = "#"
1596          ierrn = 1 ! not zero
1597!
1598!--       Search for the next line consisting of "# time",
1599!--       from there onwards the profiles will be read
1600          DO WHILE ( .NOT. ( hash == "#" .AND. ierrn == 0 ) ) 
1601         
1602            READ ( finput, *, IOSTAT=ierrn ) hash, timenudge(nt)
1603            IF ( ierrn < 0 )  EXIT rloop
1604
1605          ENDDO
1606
1607          ierrn = 0
1608          READ ( finput, *, IOSTAT=ierrn ) lowheight, lowtnudge, lowunudge,   &
1609                                           lowvnudge, lowwnudge , lowptnudge, &
1610                                           lowqnudge
1611
1612          IF ( ierrn /= 0 )  THEN
1613             message_string = 'errors in file NUDGING_DATA'
1614             CALL message( 'nudging', 'PA0366', 1, 2, 0, 6, 0 )
1615          ENDIF
1616
1617          ierrn = 0
1618          READ ( finput, *, IOSTAT=ierrn ) highheight, hightnudge, highunudge,   &
1619                                           highvnudge, highwnudge , highptnudge, &
1620                                           highqnudge
1621
1622          IF ( ierrn /= 0 )  THEN
1623             message_string = 'errors in file NUDGING_DATA'
1624             CALL message( 'nudging', 'PA0366', 1, 2, 0, 6, 0 )
1625          ENDIF
1626
1627          DO  k = nzb, nzt+1
1628             DO WHILE ( highheight < zu(k) )
1629                lowheight  = highheight
1630                lowtnudge  = hightnudge
1631                lowunudge  = highunudge
1632                lowvnudge  = highvnudge
1633                lowwnudge  = highwnudge
1634                lowptnudge = highptnudge
1635                lowqnudge  = highqnudge
1636 
1637                ierrn = 0
1638                READ ( finput, *, IOSTAT=ierrn )  highheight , hightnudge ,    &
1639                                                  highunudge , highvnudge ,    &
1640                                                  highwnudge , highptnudge,    &
1641                                                  highqnudge
1642                IF (ierrn /= 0 )  THEN
1643                   WRITE( message_string, * ) 'zu(',k,') = ', zu(k), 'm is ',  &
1644                        'higher than the maximum height in NUDING_DATA which ',&
1645                        'is ', lowheight, 'm. Interpolation on PALM ',         &
1646                        'grid is not possible.'
1647                   CALL message( 'nudging', 'PA0364', 1, 2, 0, 6, 0 )
1648                ENDIF
1649             ENDDO
1650
1651!
1652!--          Interpolation of prescribed profiles in space
1653
1654             fac = ( highheight - zu(k) ) / ( highheight - lowheight )
1655
1656             tnudge(k,nt)  = fac * lowtnudge  + ( 1.0_wp - fac ) * hightnudge
1657             unudge(k,nt)  = fac * lowunudge  + ( 1.0_wp - fac ) * highunudge
1658             vnudge(k,nt)  = fac * lowvnudge  + ( 1.0_wp - fac ) * highvnudge
1659             wnudge(k,nt)  = fac * lowwnudge  + ( 1.0_wp - fac ) * highwnudge
1660             ptnudge(k,nt) = fac * lowptnudge + ( 1.0_wp - fac ) * highptnudge
1661             qnudge(k,nt)  = fac * lowqnudge  + ( 1.0_wp - fac ) * highqnudge
1662          ENDDO
1663
1664       ENDDO rloop
1665
1666       CLOSE ( finput )
1667
1668!
1669!--    Overwrite initial profiles in case of nudging
1670       IF ( nudging )  THEN
1671          pt_init = ptnudge(:,1)
1672          u_init  = unudge(:,1)
1673          v_init  = vnudge(:,1)
1674          IF ( humidity  )  THEN ! is passive_scalar correct???
1675             q_init = qnudge(:,1)
1676          ENDIF
1677
1678          WRITE( message_string, * ) 'Initial profiles of u, v, pt and q ',    &
1679                                     'from NUDGING_DATA are used.'
1680          CALL message( 'large_scale_forcing_nudging', 'PA0370', 0, 0, 0, 6, 0 )
1681       ENDIF
1682
1683
1684    END SUBROUTINE nudge_init
1685
1686!------------------------------------------------------------------------------!
1687! Description:
1688! ------------
1689!> @todo Missing subroutine description.
1690!------------------------------------------------------------------------------!
1691    SUBROUTINE calc_tnudge ( time )
1692
1693       IMPLICIT NONE
1694
1695
1696       REAL(wp) ::  dtm         !<
1697       REAL(wp) ::  dtp         !<
1698       REAL(wp) ::  time        !<
1699
1700       INTEGER(iwp) ::  k   !<
1701       INTEGER(iwp) ::  nt  !<
1702
1703       nt = 1
1704       DO WHILE ( time > timenudge(nt) )
1705         nt = nt+1
1706       ENDDO
1707       IF ( time /= timenudge(1) ) THEN
1708         nt = nt-1
1709       ENDIF
1710
1711       dtm = ( time - timenudge(nt) ) / ( timenudge(nt+1) - timenudge(nt) )
1712       dtp = ( timenudge(nt+1) - time ) / ( timenudge(nt+1) - timenudge(nt) )
1713
1714       DO  k = nzb, nzt
1715          tmp_tnudge(k) = MAX( dt_3d, tnudge(k,nt) * dtp + tnudge(k,nt+1) * dtm )
1716       ENDDO
1717
1718    END SUBROUTINE calc_tnudge
1719
1720!------------------------------------------------------------------------------!
1721! Description:
1722! ------------
1723!> Call for all grid points
1724!------------------------------------------------------------------------------!
1725    SUBROUTINE nudge ( time, prog_var )
1726
1727       IMPLICIT NONE
1728
1729       CHARACTER (LEN=*) ::  prog_var  !<
1730
1731       REAL(wp) ::  tmp_tend    !<
1732       REAL(wp) ::  dtm         !<
1733       REAL(wp) ::  dtp         !<
1734       REAL(wp) ::  time        !<
1735
1736       INTEGER(iwp) ::  i  !<
1737       INTEGER(iwp) ::  j  !<
1738       INTEGER(iwp) ::  k  !<
1739       INTEGER(iwp) ::  nt  !<
1740
1741
1742       nt = 1
1743       DO WHILE ( time > timenudge(nt) )
1744         nt = nt+1
1745       ENDDO
1746       IF ( time /= timenudge(1) ) THEN
1747         nt = nt-1
1748       ENDIF
1749
1750       dtm = ( time - timenudge(nt) ) / ( timenudge(nt+1) - timenudge(nt) )
1751       dtp = ( timenudge(nt+1) - time ) / ( timenudge(nt+1) - timenudge(nt) )
1752
1753       SELECT CASE ( prog_var )
1754
1755          CASE ( 'u' )
1756
1757             DO  i = nxl, nxr
1758                DO  j = nys, nyn
1759
1760                   DO  k = nzb+1, nzt
1761
1762                      tmp_tend = - ( hom(k,1,1,0) - ( unudge(k,nt) * dtp +     &
1763                                     unudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
1764
1765                      tend(k,j,i) = tend(k,j,i) + tmp_tend *                   &
1766                                        MERGE( 1.0_wp, 0.0_wp,                 &
1767                                               BTEST( wall_flags_0(k,j,i), 1 ) )
1768
1769                      sums_ls_l(k,6) = sums_ls_l(k,6) + tmp_tend *             &
1770                                     weight_substep(intermediate_timestep_count)
1771                   ENDDO
1772                 
1773                   sums_ls_l(nzt+1,6) = sums_ls_l(nzt,6)
1774 
1775                ENDDO
1776            ENDDO
1777
1778          CASE ( 'v' )
1779
1780             DO  i = nxl, nxr
1781                DO  j = nys, nyn
1782
1783                   DO  k = nzb+1, nzt
1784
1785                      tmp_tend = - ( hom(k,1,2,0) - ( vnudge(k,nt) * dtp +     &
1786                                     vnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
1787
1788                      tend(k,j,i) = tend(k,j,i) + tmp_tend *                   &
1789                                        MERGE( 1.0_wp, 0.0_wp,                 &
1790                                               BTEST( wall_flags_0(k,j,i), 2 ) )
1791
1792                      sums_ls_l(k,7) = sums_ls_l(k,7) + tmp_tend *             &
1793                                     weight_substep(intermediate_timestep_count)
1794                   ENDDO
1795                 
1796                   sums_ls_l(nzt+1,7) = sums_ls_l(nzt,7)
1797
1798                ENDDO
1799            ENDDO
1800
1801          CASE ( 'pt' )
1802
1803             DO  i = nxl, nxr
1804                DO  j = nys, nyn
1805
1806                   DO  k = nzb+1, nzt
1807
1808                      tmp_tend = - ( hom(k,1,4,0) - ( ptnudge(k,nt) * dtp +    &
1809                                     ptnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
1810
1811                      tend(k,j,i) = tend(k,j,i) + tmp_tend *                   &
1812                                        MERGE( 1.0_wp, 0.0_wp,                 &
1813                                               BTEST( wall_flags_0(k,j,i), 0 ) )
1814
1815                      sums_ls_l(k,4) = sums_ls_l(k,4) + tmp_tend *             &
1816                                     weight_substep(intermediate_timestep_count)
1817                   ENDDO
1818
1819                   sums_ls_l(nzt+1,4) = sums_ls_l(nzt,4)
1820
1821                ENDDO
1822            ENDDO
1823
1824          CASE ( 'q' )
1825
1826             DO  i = nxl, nxr
1827                DO  j = nys, nyn
1828
1829                   DO  k = nzb+1, nzt
1830
1831                      tmp_tend = - ( hom(k,1,41,0) - ( qnudge(k,nt) * dtp +    &
1832                                     qnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
1833
1834                      tend(k,j,i) = tend(k,j,i) + tmp_tend *                   &
1835                                        MERGE( 1.0_wp, 0.0_wp,                 &
1836                                               BTEST( wall_flags_0(k,j,i), 0 ) )
1837
1838                      sums_ls_l(k,5) = sums_ls_l(k,5) + tmp_tend *             &
1839                                     weight_substep(intermediate_timestep_count)
1840                   ENDDO
1841                 
1842                   sums_ls_l(nzt+1,5) = sums_ls_l(nzt,5)
1843
1844                ENDDO
1845            ENDDO
1846
1847          CASE DEFAULT
1848             message_string = 'unknown prognostic variable "' // prog_var // '"'
1849             CALL message( 'nudge', 'PA0367', 1, 2, 0, 6, 0 )
1850
1851       END SELECT
1852
1853    END SUBROUTINE nudge
1854
1855
1856!------------------------------------------------------------------------------!
1857! Description:
1858! ------------
1859!> Call for grid point i,j
1860!------------------------------------------------------------------------------!
1861
1862    SUBROUTINE nudge_ij( i, j, time, prog_var )
1863
1864       IMPLICIT NONE
1865
1866
1867       CHARACTER (LEN=*) ::  prog_var  !<
1868
1869       REAL(wp) ::  tmp_tend    !<
1870       REAL(wp) ::  dtm         !<
1871       REAL(wp) ::  dtp         !<
1872       REAL(wp) ::  time        !<
1873
1874       INTEGER(iwp) ::  i  !<
1875       INTEGER(iwp) ::  j  !<
1876       INTEGER(iwp) ::  k  !<
1877       INTEGER(iwp) ::  nt  !<
1878
1879
1880       nt = 1
1881       DO WHILE ( time > timenudge(nt) )
1882         nt = nt+1
1883       ENDDO
1884       IF ( time /= timenudge(1) )  THEN
1885         nt = nt-1
1886       ENDIF
1887
1888       dtm = ( time - timenudge(nt) ) / ( timenudge(nt+1) - timenudge(nt) )
1889       dtp = ( timenudge(nt+1) - time ) / ( timenudge(nt+1) - timenudge(nt) )
1890
1891       SELECT CASE ( prog_var )
1892
1893          CASE ( 'u' )
1894
1895             DO  k = nzb+1, nzt
1896
1897                tmp_tend = - ( hom(k,1,1,0) - ( unudge(k,nt) * dtp +           &
1898                               unudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
1899
1900                tend(k,j,i) = tend(k,j,i) + tmp_tend *                         &
1901                                        MERGE( 1.0_wp, 0.0_wp,                 &
1902                                               BTEST( wall_flags_0(k,j,i), 1 ) )
1903
1904                sums_ls_l(k,6) = sums_ls_l(k,6) + tmp_tend                     &
1905                                 * weight_substep(intermediate_timestep_count)
1906             ENDDO
1907
1908             sums_ls_l(nzt+1,6) = sums_ls_l(nzt,6)
1909
1910          CASE ( 'v' )
1911
1912             DO  k = nzb+1, nzt
1913
1914                tmp_tend = - ( hom(k,1,2,0) - ( vnudge(k,nt) * dtp +           &
1915                               vnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
1916
1917                tend(k,j,i) = tend(k,j,i) + tmp_tend *                         &
1918                                        MERGE( 1.0_wp, 0.0_wp,                 &
1919                                               BTEST( wall_flags_0(k,j,i), 2 ) )
1920
1921                sums_ls_l(k,7) = sums_ls_l(k,7) + tmp_tend                     &
1922                                 * weight_substep(intermediate_timestep_count)
1923             ENDDO
1924
1925             sums_ls_l(nzt+1,7) = sums_ls_l(nzt,7)
1926
1927          CASE ( 'pt' )
1928
1929             DO  k = nzb+1, nzt
1930
1931                tmp_tend = - ( hom(k,1,4,0) - ( ptnudge(k,nt) * dtp +          &
1932                               ptnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
1933
1934                tend(k,j,i) = tend(k,j,i) + tmp_tend *                         &
1935                                        MERGE( 1.0_wp, 0.0_wp,                 &
1936                                               BTEST( wall_flags_0(k,j,i), 0 ) )
1937
1938                sums_ls_l(k,4) = sums_ls_l(k,4) + tmp_tend                     &
1939                                 * weight_substep(intermediate_timestep_count)
1940             ENDDO
1941
1942             sums_ls_l(nzt+1,4) = sums_ls_l(nzt,4)
1943
1944
1945          CASE ( 'q' )
1946
1947             DO  k = nzb+1, nzt
1948
1949                tmp_tend = - ( hom(k,1,41,0) - ( qnudge(k,nt) * dtp +          &
1950                               qnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
1951
1952                tend(k,j,i) = tend(k,j,i) + tmp_tend *                         &
1953                                        MERGE( 1.0_wp, 0.0_wp,                 &
1954                                               BTEST( wall_flags_0(k,j,i), 0 ) )
1955
1956                sums_ls_l(k,5) = sums_ls_l(k,5) + tmp_tend                     &
1957                                 * weight_substep(intermediate_timestep_count)
1958             ENDDO
1959
1960             sums_ls_l(nzt+1,5) = sums_ls_l(nzt,5)
1961
1962          CASE DEFAULT
1963             message_string = 'unknown prognostic variable "' // prog_var // '"'
1964             CALL message( 'nudge', 'PA0367', 1, 2, 0, 6, 0 )
1965
1966       END SELECT
1967
1968
1969    END SUBROUTINE nudge_ij
1970
1971
1972!------------------------------------------------------------------------------!
1973! Description:
1974! ------------
1975!> @todo Missing subroutine description.
1976!------------------------------------------------------------------------------!
1977    SUBROUTINE nudge_ref ( time )
1978
1979       IMPLICIT NONE
1980
1981       INTEGER(iwp) ::  nt                    !<
1982
1983       REAL(wp)             ::  fac           !<
1984       REAL(wp), INTENT(in) ::  time          !<
1985
1986!
1987!--    Interpolation in time of NUDGING_DATA for pt_init and q_init. This is
1988!--    needed for correct upper boundary conditions for pt and q and in case that
1989!      large scale subsidence as well as scalar Rayleigh-damping are used
1990       nt = 1
1991       DO WHILE ( time > time_vert(nt) )
1992          nt = nt + 1
1993       ENDDO
1994       IF ( time /= time_vert(nt) )  THEN
1995        nt = nt - 1
1996       ENDIF
1997
1998       fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
1999
2000       pt_init = ptnudge(:,nt) + fac * ( ptnudge(:,nt+1) - ptnudge(:,nt) )
2001       q_init  = qnudge(:,nt) + fac * ( qnudge(:,nt+1) - qnudge(:,nt) )
2002       u_init  = unudge(:,nt) + fac * ( unudge(:,nt+1) - unudge(:,nt) )
2003       v_init  = vnudge(:,nt) + fac * ( vnudge(:,nt+1) - vnudge(:,nt) )
2004
2005    END SUBROUTINE nudge_ref
2006
2007
2008 END MODULE lsf_nudging_mod
Note: See TracBrowser for help on using the repository browser.