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

Last change on this file since 3241 was 3241, checked in by raasch, 3 years ago

various changes to avoid compiler warnings (mainly removal of unused variables)

  • Property svn:keywords set to Id
File size: 78.1 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 3241 2018-09-12 15:02:00Z raasch $
27! unused variables removed
28!
29! 3183 2018-07-27 14:25:55Z suehring
30! * Adjustment to new Inifor version:
31!   - No vertical interpolation/extrapolation of lateral boundary data required
32!     any more (Inifor can treat grid stretching now
33!   - Revise initialization in case of COSMO forcing
34! * Rename variables and subroutines for offline nesting
35!
36! 3182 2018-07-27 13:36:03Z suehring
37! Error messages revised
38!
39! 3045 2018-05-28 07:55:41Z Giersch
40! Error messages revised
41!
42! 3026 2018-05-22 10:30:53Z schwenkel
43! Changed the name specific humidity to mixing ratio, since we are computing
44! mixing ratios.
45!
46! 2970 2018-04-13 15:09:23Z suehring
47! Bugfix in old large-scale forcing mode
48!
49! 2938 2018-03-27 15:52:42Z suehring
50! Further improvements for nesting in larger-scale model
51!
52! 2863 2018-03-08 11:36:25Z suehring
53! Corrected "Former revisions" section
54!
55! 2696 2017-12-14 17:12:51Z kanani
56! Change in file header (GPL part)
57! Forcing with larger-scale models implemented (MS)
58!
59! 2342 2017-08-08 11:00:43Z boeske
60! fixed check if surface forcing data is available until end of simulation
61!
62! 2320 2017-07-21 12:47:43Z suehring
63! initial revision
64!
65! Description:
66! ------------
67!> Calculates large scale forcings (geostrophic wind and subsidence velocity) as
68!> well as surfaces fluxes dependent on time given in an external file (LSF_DATA).
69!> Moreover, module contains nudging routines, where u, v, pt and q are nudged
70!> to given profiles on a relaxation timescale tnudge.
71!> Profiles are read in from NUDGING_DATA.
72!> Code is based on Neggers et al. (2012) and also in parts on DALES and UCLA-LES.
73!> @todo: Revise reading of ASCII-files
74!> @todo: Remove unused variables and control flags
75!> @todo: Revise large-scale facing of surface variables
76!> @todo: Revise control flags lsf_exception, lsf_surf, lsf_vert, etc.
77!--------------------------------------------------------------------------------!
78 MODULE lsf_nudging_mod
79
80    USE arrays_3d,                                                             &
81        ONLY:  dzw, e, diss, heatflux_input_conversion, pt, pt_init, q,        &
82               q_init, s, tend, u, u_init, ug, v, v_init, vg, w, w_subs,       &
83               waterflux_input_conversion, zu, zw                 
84
85    USE control_parameters,                                                    &
86        ONLY:  bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, bc_dirichlet_s, &
87               bc_lr, bc_ns, bc_pt_b, bc_q_b, constant_diffusion,              &
88               constant_heatflux, constant_waterflux,                          &
89               data_output_pr, dt_3d, end_time,                                &
90               humidity, initializing_actions, intermediate_timestep_count,    &
91               ibc_pt_b, ibc_q_b,                                              &
92               large_scale_forcing, large_scale_subsidence, lsf_surf, lsf_vert,&
93               lsf_exception, message_string, nesting_offline, neutral,        &
94               nudging, passive_scalar, pt_surface, ocean, q_surface,          &
95               surface_heatflux, surface_pressure, surface_waterflux,          &
96               topography, use_subsidence_tendencies
97
98    USE grid_variables
99
100    USE indices,                                                               &
101        ONLY:  nbgp, ngp_sums_ls, nx, nxl, nxlg, nxlu, nxr, nxrg, ny, nys,     &
102               nysv, nysg, nyn, nyng, nzb, nz, nzt, wall_flags_0
103
104    USE kinds
105
106    USE netcdf_data_input_mod,                                                 &
107        ONLY:  nest_offl
108
109    USE pegrid
110
111    USE surface_mod,                                                           &
112        ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
113
114    USE statistics,                                                            &
115        ONLY:  hom, statistic_regions, sums_ls_l, weight_substep
116
117    INTEGER(iwp) ::  nlsf = 1000                       !< maximum number of profiles in LSF_DATA (large scale forcing)
118    INTEGER(iwp) ::  ntnudge = 1000                    !< maximum number of profiles in NUDGING_DATA (nudging)
119
120    REAL(wp) ::  d_area_t
121
122    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ptnudge     !< vertical profile of pot. temperature interpolated to vertical grid (nudging)
123    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  qnudge      !< vertical profile of water vapor mixing ratio interpolated to vertical grid (nudging)
124    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tnudge      !< vertical profile of nudging time scale interpolated to vertical grid (nudging) 
125    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  td_lsa_lpt  !< temperature tendency due to large scale advection (large scale forcing)
126    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  td_lsa_q    !< water vapor mixing ratio tendency due to large scale advection (large scale forcing)
127    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  td_sub_lpt  !< temperature tendency due to subsidence/ascent (large scale forcing)
128    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  td_sub_q    !< water vapor mixing ratio tendency due to subsidence/ascent (large scale forcing)
129    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ug_vert     !< vertical profile of geostrophic wind component in x-direction interpolated to vertical grid (large scale forcing)
130    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  unudge      !< vertical profile of wind component in x-direction interpolated to vertical grid (nudging)
131    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  vnudge      !< vertical profile of wind component in y-direction interpolated to vertical grid (nudging)
132    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  vg_vert     !< vertical profile of geostrophic wind component in y-direction interpolated to vertical grid (large scale forcing)
133    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  wnudge      !< vertical profile of subsidence/ascent velocity interpolated to vertical grid (nudging) ???
134    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  wsubs_vert  !< vertical profile of wind component in z-direction interpolated to vertical grid (nudging) ???
135
136    REAL(wp), DIMENSION(:), ALLOCATABLE ::  shf_surf      !< time-dependent surface sensible heat flux (large scale forcing)
137    REAL(wp), DIMENSION(:), ALLOCATABLE ::  timenudge     !< times at which vertical profiles are defined in NUDGING_DATA (nudging)
138    REAL(wp), DIMENSION(:), ALLOCATABLE ::  time_surf     !< times at which surface values/fluxes are defined in LSF_DATA (large scale forcing)
139    REAL(wp), DIMENSION(:), ALLOCATABLE ::  time_vert     !< times at which vertical profiles are defined in LSF_DATA (large scale forcing)
140
141    REAL(wp), DIMENSION(:), ALLOCATABLE ::  tmp_tnudge    !< current nudging time scale
142
143    REAL(wp), DIMENSION(:), ALLOCATABLE ::  p_surf        !< time-dependent surface pressure (large scale forcing)
144    REAL(wp), DIMENSION(:), ALLOCATABLE ::  pt_surf       !< time-dependent surface temperature (large scale forcing)
145    REAL(wp), DIMENSION(:), ALLOCATABLE ::  qsws_surf     !< time-dependent surface latent heat flux (large scale forcing)
146    REAL(wp), DIMENSION(:), ALLOCATABLE ::  q_surf        !< time-dependent surface water vapor mixing ratio (large scale forcing)
147
148    SAVE
149    PRIVATE
150!
151!-- Public subroutines
152    PUBLIC calc_tnudge, ls_forcing_surf, ls_forcing_vert, ls_advec, lsf_init,  &
153           lsf_nudging_check_parameters, nudge_init,                           &
154           lsf_nudging_check_data_output_pr, lsf_nudging_header,               &
155           lsf_nesting_offline, lsf_nesting_offline_mass_conservation,         & 
156           nudge, nudge_ref
157           
158!
159!-- Public variables
160    PUBLIC qsws_surf, shf_surf, td_lsa_lpt, td_lsa_q, td_sub_lpt,              &
161           td_sub_q, time_vert
162
163
164    INTERFACE ls_advec
165       MODULE PROCEDURE ls_advec
166       MODULE PROCEDURE ls_advec_ij
167    END INTERFACE ls_advec
168
169    INTERFACE nudge
170       MODULE PROCEDURE nudge
171       MODULE PROCEDURE nudge_ij
172    END INTERFACE nudge
173
174 CONTAINS
175
176
177!------------------------------------------------------------------------------!
178! Description:
179! ------------
180!> In this subroutine a constant mass within the model domain is guaranteed.
181!> Larger-scale models may be based on a compressible equation system, which is
182!> not consistent with PALMs incompressible equation system. In order to avoid
183!> a decrease or increase of mass during the simulation, non-divergent flow
184!> through the lateral and top boundaries is compensated by the vertical wind
185!> component at the top boundary.
186!------------------------------------------------------------------------------!
187    SUBROUTINE lsf_nesting_offline_mass_conservation
188
189       USE control_parameters,                                                 &
190           ONLY:  volume_flow
191
192       IMPLICIT NONE
193
194       INTEGER(iwp) ::  i !< grid index in x-direction
195       INTEGER(iwp) ::  j !< grid index in y-direction
196       INTEGER(iwp) ::  k !< grid index in z-direction
197
198       REAL(wp) ::  w_correct                       !< vertical velocity increment required to compensate non-divergent flow through the boundaries
199       REAL(wp), DIMENSION(1:3) ::  volume_flow_l   !< local volume flow
200
201       volume_flow   = 0.0_wp
202       volume_flow_l = 0.0_wp
203
204       d_area_t = 1.0_wp / ( ( nx + 1 ) * dx * ( ny + 1 ) * dy )
205
206       IF ( bc_dirichlet_l )  THEN
207          i = nxl
208          DO  j = nys, nyn
209             DO  k = nzb+1, nzt
210                volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k) * dy   &
211                                     * MERGE( 1.0_wp, 0.0_wp,                  &
212                                              BTEST( wall_flags_0(k,j,i), 1 ) )
213             ENDDO
214          ENDDO
215       ENDIF
216       IF ( bc_dirichlet_r )  THEN
217          i = nxr+1
218          DO  j = nys, nyn
219             DO  k = nzb+1, nzt
220                volume_flow_l(1) = volume_flow_l(1) - u(k,j,i) * dzw(k) * dy   &
221                                     * MERGE( 1.0_wp, 0.0_wp,                  &
222                                              BTEST( wall_flags_0(k,j,i), 1 ) )
223             ENDDO
224          ENDDO
225       ENDIF
226       IF ( bc_dirichlet_s )  THEN
227          j = nys
228          DO  i = nxl, nxr
229             DO  k = nzb+1, nzt
230                volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k) * dx   &
231                                     * MERGE( 1.0_wp, 0.0_wp,                  &
232                                              BTEST( wall_flags_0(k,j,i), 2 ) )
233             ENDDO
234          ENDDO
235       ENDIF
236       IF ( bc_dirichlet_n )  THEN
237          j = nyn+1
238          DO  i = nxl, nxr
239             DO  k = nzb+1, nzt
240                volume_flow_l(2) = volume_flow_l(2) - v(k,j,i) * dzw(k) * dx   &
241                                     * MERGE( 1.0_wp, 0.0_wp,                  &
242                                              BTEST( wall_flags_0(k,j,i), 2 ) )
243             ENDDO
244          ENDDO
245       ENDIF
246!
247!--    Top boundary
248       k = nzt
249       DO  i = nxl, nxr
250          DO  j = nys, nyn
251             volume_flow_l(3) = volume_flow_l(3) - w(k,j,i) * dx * dy
252          ENDDO
253       ENDDO
254
255#if defined( __parallel )
256       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
257       CALL MPI_ALLREDUCE( volume_flow_l, volume_flow, 3, MPI_REAL, MPI_SUM,      &
258                           comm2d, ierr )
259#else
260       volume_flow = volume_flow_l
261#endif
262
263       w_correct = SUM( volume_flow ) * d_area_t
264
265       DO  i = nxl, nxr
266          DO  j = nys, nyn
267             DO  k = nzt, nzt + 1
268                w(k,j,i) = w(k,j,i) + w_correct
269             ENDDO
270          ENDDO
271       ENDDO
272
273    END SUBROUTINE lsf_nesting_offline_mass_conservation
274
275
276!------------------------------------------------------------------------------!
277! Description:
278! ------------
279!> Set the lateral and top boundary conditions in case the PALM domain is
280!> nested offline in a mesoscale model.
281!------------------------------------------------------------------------------!
282    SUBROUTINE lsf_nesting_offline
283
284       USE control_parameters,                                                 &
285           ONLY:  bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r,              &
286                  bc_dirichlet_s, humidity, neutral, passive_scalar, rans_mode,&
287                  rans_tke_e, time_since_reference_point                     
288
289       IMPLICIT NONE
290
291       INTEGER(iwp) ::  i !< running index x-direction
292       INTEGER(iwp) ::  j !< running index y-direction
293       INTEGER(iwp) ::  k !< running index z-direction
294
295       REAL(wp) ::  ddt_lsf !< inverse value of time resolution of forcing data
296       REAL(wp) ::  t_ref   !< time past since last reference step
297     
298!
299!--    Calculate time interval of forcing data       
300       ddt_lsf = 1.0_wp / ( nest_offl%time(nest_offl%tind_p) -                 &
301                            nest_offl%time(nest_offl%tind) )
302!
303!--    Calculate reziproke time past since last reference step. Please note,
304!--    the time coordinate is still not updated, so that the actual time need
305!--    to be incremented by dt_3d. Moreover, note that the simulation time
306!--    passed since simulation start is time_since_reference_point, not
307!--    simulated_time!
308       t_ref = time_since_reference_point + dt_3d -                            &
309                                            nest_offl%time(nest_offl%tind)
310                                           
311       IF ( bc_dirichlet_l )  THEN
312
313          DO  j = nys, nyn
314             DO  k = nzb+1, nzt
315                u(k,j,nxlg:nxl) = nest_offl%u_left(0,k,j) + ddt_lsf * t_ref  * &
316                       ( nest_offl%u_left(1,k,j) - nest_offl%u_left(0,k,j) ) * &
317                         MERGE( 1.0_wp, 0.0_wp,                                &
318                                BTEST( wall_flags_0(k,j,nxlg:nxl), 1 ) )
319             ENDDO
320          ENDDO
321
322          DO  j = nys, nyn
323             DO  k = nzb+1, nzt-1
324                w(k,j,nxlg:nxl-1) = nest_offl%w_left(0,k,j) + ddt_lsf * t_ref *&
325                       ( nest_offl%w_left(1,k,j) - nest_offl%w_left(0,k,j) )  *&
326                         MERGE( 1.0_wp, 0.0_wp,                                &
327                                BTEST( wall_flags_0(k,j,nxlg:nxl-1), 3 ) )
328             ENDDO
329          ENDDO
330
331          DO  j = nysv, nyn
332             DO  k = nzb+1, nzt
333                v(k,j,nxlg:nxl-1) = nest_offl%v_left(0,k,j) + ddt_lsf * t_ref *&
334                       ( nest_offl%v_left(1,k,j) - nest_offl%v_left(0,k,j) )  *&
335                         MERGE( 1.0_wp, 0.0_wp,                                &
336                                BTEST( wall_flags_0(k,j,nxlg:nxl-1), 2 ) )
337             ENDDO
338          ENDDO
339
340          IF ( .NOT. neutral )  THEN
341             DO  j = nys, nyn
342                DO  k = nzb+1, nzt
343                   pt(k,j,nxlg:nxl-1) = nest_offl%pt_left(0,k,j) + ddt_lsf *   &
344                                                                   t_ref   *   &
345                       ( nest_offl%pt_left(1,k,j) - nest_offl%pt_left(0,k,j) )
346 
347                ENDDO
348             ENDDO
349          ENDIF
350
351          IF ( humidity )  THEN
352             DO  j = nys, nyn
353                DO  k = nzb+1, nzt
354                   q(k,j,nxlg:nxl-1) = nest_offl%q_left(0,k,j) + ddt_lsf   *   &
355                                                                 t_ref     *   &
356                       ( nest_offl%q_left(1,k,j) - nest_offl%q_left(0,k,j) )
357 
358                ENDDO
359             ENDDO
360          ENDIF
361
362       ENDIF
363
364       IF ( bc_dirichlet_r )  THEN
365
366          DO  j = nys, nyn
367             DO  k = nzb+1, nzt
368                u(k,j,nxr+1:nxrg) = nest_offl%u_right(0,k,j) + ddt_lsf * t_ref *&
369                      ( nest_offl%u_right(1,k,j) - nest_offl%u_right(0,k,j) )  *&
370                         MERGE( 1.0_wp, 0.0_wp,                                 &
371                                BTEST( wall_flags_0(k,j,nxr+1:nxrg), 1 ) )
372
373             ENDDO
374          ENDDO
375          DO  j = nys, nyn
376             DO  k = nzb+1, nzt-1
377                w(k,j,nxr+1:nxrg) = nest_offl%w_right(0,k,j) + ddt_lsf * t_ref *&
378                      ( nest_offl%w_right(1,k,j) - nest_offl%w_right(0,k,j) )  *&
379                         MERGE( 1.0_wp, 0.0_wp,                                 &
380                                BTEST( wall_flags_0(k,j,nxr+1:nxrg), 3 ) )
381             ENDDO
382          ENDDO
383
384          DO  j = nysv, nyn
385             DO  k = nzb+1, nzt
386                v(k,j,nxr+1:nxrg) = nest_offl%v_right(0,k,j) + ddt_lsf * t_ref *&
387                      ( nest_offl%v_right(1,k,j) - nest_offl%v_right(0,k,j) )  *&
388                         MERGE( 1.0_wp, 0.0_wp,                                 &
389                                BTEST( wall_flags_0(k,j,nxr+1:nxrg), 2 ) )
390             ENDDO
391          ENDDO
392
393          IF ( .NOT. neutral )  THEN
394             DO  j = nys, nyn
395                DO  k = nzb+1, nzt
396                   pt(k,j,nxr+1:nxrg) = nest_offl%pt_right(0,k,j) + ddt_lsf *  &
397                                                                    t_ref   *  &
398                     ( nest_offl%pt_right(1,k,j) - nest_offl%pt_right(0,k,j) )
399 
400                ENDDO
401             ENDDO
402          ENDIF
403
404          IF ( humidity )  THEN
405             DO  j = nys, nyn
406                DO  k = nzb+1, nzt
407                   q(k,j,nxr+1:nxrg) = nest_offl%q_right(0,k,j) + ddt_lsf   *  &
408                                                                    t_ref   *  &
409                       ( nest_offl%q_right(1,k,j) - nest_offl%q_right(0,k,j) )
410 
411                ENDDO
412             ENDDO
413          ENDIF
414
415       ENDIF
416
417       IF ( bc_dirichlet_s )  THEN
418
419          DO  i = nxl, nxr
420             DO  k = nzb+1, nzt
421                v(k,nysg:nys,i)   = nest_offl%v_south(0,k,i) + ddt_lsf * t_ref *&
422                      ( nest_offl%v_south(1,k,i) - nest_offl%v_south(0,k,i) )  *&
423                         MERGE( 1.0_wp, 0.0_wp,                                 &
424                                BTEST( wall_flags_0(k,nysg:nys,i), 2 ) )
425             ENDDO
426          ENDDO
427
428          DO  i = nxl, nxr
429             DO  k = nzb+1, nzt-1
430                w(k,nysg:nys-1,i) = nest_offl%w_south(0,k,i) + ddt_lsf * t_ref  *&
431                        ( nest_offl%w_south(1,k,i) - nest_offl%w_south(0,k,i) ) *&
432                           MERGE( 1.0_wp, 0.0_wp,                                &
433                                  BTEST( wall_flags_0(k,nysg:nys-1,i), 3 ) )
434             ENDDO
435          ENDDO
436
437          DO  i = nxlu, nxr
438             DO  k = nzb+1, nzt
439                u(k,nysg:nys-1,i) = nest_offl%u_south(0,k,i) + ddt_lsf * t_ref  *&
440                        ( nest_offl%u_south(1,k,i) - nest_offl%u_south(0,k,i) ) *&
441                           MERGE( 1.0_wp, 0.0_wp,                                &
442                                  BTEST( wall_flags_0(k,nysg:nys-1,i), 1 ) )
443             ENDDO
444          ENDDO
445
446          IF ( .NOT. neutral )  THEN
447             DO  i = nxl, nxr
448                DO  k = nzb+1, nzt
449                   pt(k,nysg:nys-1,i) = nest_offl%pt_south(0,k,i) + ddt_lsf *  &
450                                                                    t_ref   *  &
451                     ( nest_offl%pt_south(1,k,i) - nest_offl%pt_south(0,k,i) )
452 
453                ENDDO
454             ENDDO
455          ENDIF
456
457          IF ( humidity )  THEN
458             DO  i = nxl, nxr
459                DO  k = nzb+1, nzt
460                   q(k,nysg:nys-1,i) = nest_offl%q_south(0,k,i) + ddt_lsf   *  &
461                                                                    t_ref   *  &
462                       ( nest_offl%q_south(1,k,i) - nest_offl%q_south(0,k,i) )
463 
464                ENDDO
465             ENDDO
466          ENDIF
467
468       ENDIF
469
470       IF ( bc_dirichlet_n )  THEN
471
472          DO  i = nxl, nxr
473             DO  k = nzb+1, nzt
474                v(k,nyn+1:nyng,i)   = nest_offl%v_north(0,k,i) + ddt_lsf * t_ref *&
475                        ( nest_offl%v_north(1,k,i) - nest_offl%v_north(0,k,i) )  *&
476                           MERGE( 1.0_wp, 0.0_wp,                                 &
477                                  BTEST( wall_flags_0(k,nyn+1:nyng,i), 2 ) )
478             ENDDO
479          ENDDO
480          DO  i = nxl, nxr
481             DO  k = nzb+1, nzt-1
482                w(k,nyn+1:nyng,i) = nest_offl%w_north(0,k,i) + ddt_lsf * t_ref  *&
483                        ( nest_offl%w_north(1,k,i) - nest_offl%w_north(0,k,i) ) *&
484                           MERGE( 1.0_wp, 0.0_wp,                                &
485                                  BTEST( wall_flags_0(k,nyn+1:nyng,i), 3 ) )
486             ENDDO
487          ENDDO
488
489          DO  i = nxlu, nxr
490             DO  k = nzb+1, nzt
491                u(k,nyn+1:nyng,i) = nest_offl%u_north(0,k,i) + ddt_lsf * t_ref  *&
492                        ( nest_offl%u_north(1,k,i) - nest_offl%u_north(0,k,i) ) *&
493                           MERGE( 1.0_wp, 0.0_wp,                                &
494                                  BTEST( wall_flags_0(k,nyn+1:nyng,i), 1 ) )
495
496             ENDDO
497          ENDDO
498
499          IF ( .NOT. neutral )  THEN
500             DO  i = nxl, nxr
501                DO  k = nzb+1, nzt
502                   pt(k,nyn+1:nyng,i) = nest_offl%pt_north(0,k,i) + ddt_lsf *  &
503                                                                    t_ref   *  &
504                     ( nest_offl%pt_north(1,k,i) - nest_offl%pt_north(0,k,i) )
505 
506                ENDDO
507             ENDDO
508          ENDIF
509
510          IF ( humidity )  THEN
511             DO  i = nxl, nxr
512                DO  k = nzb+1, nzt
513                   q(k,nyn+1:nyng,i) = nest_offl%q_north(0,k,i) + ddt_lsf   *  &
514                                                                    t_ref   *  &
515                       ( nest_offl%q_north(1,k,i) - nest_offl%q_north(0,k,i) )
516 
517                ENDDO
518             ENDDO
519          ENDIF
520
521       ENDIF
522!
523!--    Top boundary.
524       DO  i = nxlu, nxr
525          DO  j = nys, nyn
526             u(nzt+1,j,i) = nest_offl%u_top(0,j,i) + ddt_lsf * t_ref *         &
527                        ( nest_offl%u_top(1,j,i) - nest_offl%u_top(0,j,i) ) *  &
528                           MERGE( 1.0_wp, 0.0_wp,                              &
529                                  BTEST( wall_flags_0(nzt+1,j,i), 1 ) )
530          ENDDO
531       ENDDO
532
533       DO  i = nxl, nxr
534          DO  j = nysv, nyn
535             v(nzt+1,j,i) = nest_offl%v_top(0,j,i) + ddt_lsf * t_ref *         &
536                        ( nest_offl%v_top(1,j,i) - nest_offl%v_top(0,j,i) ) *  &
537                           MERGE( 1.0_wp, 0.0_wp,                              &
538                                  BTEST( wall_flags_0(nzt+1,j,i), 2 ) )
539          ENDDO
540       ENDDO
541
542       DO  i = nxl, nxr
543          DO  j = nys, nyn
544             w(nzt:nzt+1,j,i) = nest_offl%w_top(0,j,i) + ddt_lsf * t_ref *     &
545                        ( nest_offl%w_top(1,j,i) - nest_offl%w_top(0,j,i) ) *  &
546                           MERGE( 1.0_wp, 0.0_wp,                              &
547                                  BTEST( wall_flags_0(nzt:nzt+1,j,i), 3 ) )
548          ENDDO
549       ENDDO
550
551
552       IF ( .NOT. neutral )  THEN
553          DO  i = nxl, nxr
554             DO  j = nys, nyn
555                pt(nzt+1,j,i) = nest_offl%pt_top(0,j,i) + ddt_lsf * t_ref *    &
556                        ( nest_offl%pt_top(1,j,i) - nest_offl%pt_top(0,j,i) )
557             ENDDO
558          ENDDO
559       ENDIF
560
561       IF ( humidity )  THEN
562          DO  i = nxl, nxr
563             DO  j = nys, nyn
564                q(nzt+1,j,i) = nest_offl%q_top(0,j,i) + ddt_lsf * t_ref *      &
565                        ( nest_offl%q_top(1,j,i) - nest_offl%q_top(0,j,i) )
566             ENDDO
567          ENDDO
568       ENDIF
569!
570!--    At the edges( left-south, left-north, right-south and right-north) set
571!--    data on ghost points.
572       IF ( bc_dirichlet_l  .AND.  bc_dirichlet_s )  THEN
573          DO  i = 1, nbgp
574             u(:,nys-i,nxlg:nxl)   = u(:,nys,nxlg:nxl)
575             w(:,nys-i,nxlg:nxl-1) = w(:,nys,nxlg:nxl-1)
576             IF ( .NOT. neutral )  pt(:,nys-i,nxlg:nxl-1) = pt(:,nys,nxlg:nxl-1)
577             IF ( humidity      )  q(:,nys-i,nxlg:nxl-1)  = q(:,nys,nxlg:nxl-1)
578          ENDDO
579          DO  i = 1, nbgp+1
580             v(:,nysv-i,nxlg:nxl-1) = v(:,nysv,nxlg:nxl-1)
581          ENDDO
582       ENDIF
583       IF ( bc_dirichlet_l  .AND.  bc_dirichlet_n )  THEN
584          DO  i = 1, nbgp
585             u(:,nyn+i,nxlg:nxl)   = u(:,nyn,nxlg:nxl)
586             v(:,nyn+i,nxlg:nxl-1) = v(:,nyn,nxlg:nxl-1)
587             w(:,nyn+i,nxlg:nxl-1) = w(:,nyn,nxlg:nxl-1)
588             IF ( .NOT. neutral )  pt(:,nyn+i,nxlg:nxl-1) = pt(:,nyn,nxlg:nxl-1)
589             IF ( humidity      )  q(:,nyn+i,nxlg:nxl-1)  = q(:,nyn,nxlg:nxl-1)
590          ENDDO
591       ENDIF
592       IF ( bc_dirichlet_r  .AND.  bc_dirichlet_s )  THEN
593          DO  i = 1, nbgp
594             u(:,nys-i,nxr+1:nxrg) = u(:,nys,nxr+1:nxrg)
595             w(:,nys-i,nxr+1:nxrg) = w(:,nys,nxr+1:nxrg)
596             IF ( .NOT. neutral )  pt(:,nys-i,nxr+1:nxrg) = pt(:,nys,nxr+1:nxrg)
597             IF ( humidity      )  q(:,nys-i,nxr+1:nxrg)  = q(:,nys,nxr+1:nxrg)
598          ENDDO
599          DO  i = 1, nbgp+1
600             v(:,nysv-i,nxr+1:nxrg) = v(:,nysv,nxr+1:nxrg)
601          ENDDO
602       ENDIF
603       IF ( bc_dirichlet_r  .AND.  bc_dirichlet_n )  THEN
604          DO  i = 1, nbgp
605             u(:,nyn+i,nxr+1:nxrg) = u(:,nyn,nxr+1:nxrg)
606             v(:,nyn+i,nxr+1:nxrg) = v(:,nyn,nxr+1:nxrg)
607             w(:,nyn+i,nxr+1:nxrg) = w(:,nyn,nxr+1:nxrg)
608             IF ( .NOT. neutral )  pt(:,nyn+i,nxr+1:nxrg) = pt(:,nyn,nxr+1:nxrg)
609             IF ( humidity      )  q(:,nyn+i,nxr+1:nxrg)  = q(:,nyn,nxr+1:nxrg)
610          ENDDO
611       ENDIF
612!
613!--    Moreover, set Neumann boundary condition for subgrid-scale TKE,
614!--    passive scalar, dissipation, and chemical species if required
615       IF ( rans_mode  .AND.  rans_tke_e )  THEN
616          IF (  bc_dirichlet_l )  diss(:,:,nxl-1) = diss(:,:,nxl)
617          IF (  bc_dirichlet_r )  diss(:,:,nxr+1) = diss(:,:,nxr)
618          IF (  bc_dirichlet_s )  diss(:,nys-1,:) = diss(:,nys,:)
619          IF (  bc_dirichlet_n )  diss(:,nyn+1,:) = diss(:,nyn,:)
620       ENDIF
621       IF ( .NOT. constant_diffusion )  THEN
622          IF (  bc_dirichlet_l )  e(:,:,nxl-1) = e(:,:,nxl)
623          IF (  bc_dirichlet_r )  e(:,:,nxr+1) = e(:,:,nxr)
624          IF (  bc_dirichlet_s )  e(:,nys-1,:) = e(:,nys,:)
625          IF (  bc_dirichlet_n )  e(:,nyn+1,:) = e(:,nyn,:)
626          e(nzt+1,:,:) = e(nzt,:,:)
627       ENDIF
628       IF ( passive_scalar )  THEN
629          IF (  bc_dirichlet_l )  s(:,:,nxl-1) = s(:,:,nxl)
630          IF (  bc_dirichlet_r )  s(:,:,nxr+1) = s(:,:,nxr)
631          IF (  bc_dirichlet_s )  s(:,nys-1,:) = s(:,nys,:)
632          IF (  bc_dirichlet_n )  s(:,nyn+1,:) = s(:,nyn,:)
633       ENDIF
634
635
636
637       CALL exchange_horiz( u, nbgp )
638       CALL exchange_horiz( v, nbgp )
639       CALL exchange_horiz( w, nbgp )
640       IF ( .NOT. neutral )  CALL exchange_horiz( pt, nbgp )
641       IF ( humidity      )  CALL exchange_horiz( q,  nbgp )
642
643!
644!--    Set surface pressure. Please note, time-dependent surface
645!--    pressure would require changes in anelastic approximation and
646!--    treatment of fluxes.
647!--    For the moment, comment this out!
648!      surface_pressure = nest_offl%surface_pressure(nest_offl%tind) +                 &
649!                                                      ddt_lsf * t_ref *       &
650!                                    ( nest_offl%surface_pressure(nest_offl%tind_p)    &
651!                                    - nest_offl%surface_pressure(nest_offl%tind) )
652
653    END SUBROUTINE lsf_nesting_offline
654
655!------------------------------------------------------------------------------!
656! Description:
657! ------------
658!> @todo Missing subroutine description.
659!------------------------------------------------------------------------------!
660    SUBROUTINE lsf_nudging_check_parameters
661
662       IMPLICIT NONE
663!
664!--    Check nudging and large scale forcing from external file
665       IF ( nudging  .AND.  (  .NOT.  large_scale_forcing ) )  THEN
666          message_string = 'Nudging requires large_scale_forcing = .T.. &'//   &
667                        'Surface fluxes and geostrophic wind should be &'//    &
668                        'prescribed in file LSF_DATA'
669          CALL message( 'check_parameters', 'PA0374', 1, 2, 0, 6, 0 )
670       ENDIF
671
672       IF ( large_scale_forcing  .AND.  ( bc_lr /= 'cyclic'  .OR.              &
673                                          bc_ns /= 'cyclic' ) )  THEN
674          message_string = 'Non-cyclic lateral boundaries do not allow for &'//&
675                        'the usage of large scale forcing from external file.'
676          CALL message( 'check_parameters', 'PA0375', 1, 2, 0, 6, 0 )
677       ENDIF
678
679       IF ( large_scale_forcing  .AND.  (  .NOT.  humidity ) )  THEN
680          message_string = 'The usage of large scale forcing from external &'//&
681                        'file LSF_DATA requires humidity = .T..'
682          CALL message( 'check_parameters', 'PA0376', 1, 2, 0, 6, 0 )
683       ENDIF
684
685       IF ( large_scale_forcing  .AND.  passive_scalar )  THEN
686          message_string = 'The usage of large scale forcing from external &'// &
687                        'file LSF_DATA is not implemented for passive scalars'
688          CALL message( 'check_parameters', 'PA0440', 1, 2, 0, 6, 0 )
689       ENDIF
690
691       IF ( large_scale_forcing  .AND.  topography /= 'flat'                   &
692                              .AND.  .NOT.  lsf_exception )  THEN
693          message_string = 'The usage of large scale forcing from external &'//&
694                        'file LSF_DATA is not implemented for non-flat topography'
695          CALL message( 'check_parameters', 'PA0377', 1, 2, 0, 6, 0 )
696       ENDIF
697
698       IF ( large_scale_forcing  .AND.  ocean )  THEN
699          message_string = 'The usage of large scale forcing from external &'//&
700                        'file LSF_DATA is not implemented for ocean runs'
701          CALL message( 'check_parameters', 'PA0378', 1, 2, 0, 6, 0 )
702       ENDIF
703
704    END SUBROUTINE lsf_nudging_check_parameters
705
706!------------------------------------------------------------------------------!
707! Description:
708! ------------
709!> Check data output of profiles for land surface model
710!------------------------------------------------------------------------------!
711    SUBROUTINE lsf_nudging_check_data_output_pr( variable, var_count, unit,    &
712                                                 dopr_unit )
713 
714       USE profil_parameter
715
716       IMPLICIT NONE
717   
718       CHARACTER (LEN=*) ::  unit      !<
719       CHARACTER (LEN=*) ::  variable  !<
720       CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
721 
722       INTEGER(iwp) ::  var_count     !<
723
724       SELECT CASE ( TRIM( variable ) )
725       
726
727          CASE ( 'td_lsa_lpt' )
728             IF (  .NOT.  large_scale_forcing )  THEN
729                message_string = 'data_output_pr = ' //                        &
730                                 TRIM( data_output_pr(var_count) ) //          &
731                                 ' is not implemented for ' //                 &
732                                 'large_scale_forcing = .FALSE.'
733                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393',    &
734                               1, 2, 0, 6, 0 )
735             ELSE
736                dopr_index(var_count) = 81
737                dopr_unit             = 'K/s'
738                unit                  = 'K/s'
739                hom(:,2,81,:) = SPREAD( zu, 2, statistic_regions+1 )
740             ENDIF
741
742          CASE ( 'td_lsa_q' )
743             IF (  .NOT.  large_scale_forcing )  THEN
744                message_string = 'data_output_pr = ' //                        &
745                                 TRIM( data_output_pr(var_count) ) //          &
746                                 ' is not implemented for ' //                 &
747                                 'large_scale_forcing = .FALSE.'
748                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393',    &
749                               1, 2, 0, 6, 0 )
750             ELSE
751                dopr_index(var_count) = 82
752                dopr_unit             = 'kg/kgs'
753                unit                  = 'kg/kgs'
754                hom(:,2,82,:) = SPREAD( zu, 2, statistic_regions+1 )
755             ENDIF
756          CASE ( 'td_sub_lpt' )
757             IF (  .NOT.  large_scale_forcing )  THEN
758                message_string = 'data_output_pr = ' //                        &
759                                 TRIM( data_output_pr(var_count) ) //          &
760                                 ' is not implemented for ' //                 &
761                                 'large_scale_forcing = .FALSE.'
762                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393',    &
763                               1, 2, 0, 6, 0 )
764             ELSE
765                dopr_index(var_count) = 83
766                dopr_unit             = 'K/s'
767                unit                  = 'K/s'
768                hom(:,2,83,:) = SPREAD( zu, 2, statistic_regions+1 )
769             ENDIF
770
771          CASE ( 'td_sub_q' )
772             IF (  .NOT.  large_scale_forcing )  THEN
773                message_string = 'data_output_pr = ' //                        &
774                                 TRIM( data_output_pr(var_count) ) //          &
775                                 ' is not implemented for ' //                 &
776                                 'large_scale_forcing = .FALSE.'
777                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393',    &
778                               1, 2, 0, 6, 0 )
779             ELSE
780                dopr_index(var_count) = 84
781                dopr_unit             = 'kg/kgs'
782                unit                  = 'kg/kgs'
783                hom(:,2,84,:) = SPREAD( zu, 2, statistic_regions+1 )
784             ENDIF
785
786          CASE ( 'td_nud_lpt' )
787             IF (  .NOT.  nudging )  THEN
788                message_string = 'data_output_pr = ' //                        &
789                                 TRIM( data_output_pr(var_count) ) //          &
790                                 ' is not implemented for ' //                 &
791                                 'nudging = .FALSE.'
792                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0394',    &
793                               1, 2, 0, 6, 0 )
794             ELSE
795                dopr_index(var_count) = 85
796                dopr_unit             = 'K/s'
797                unit                  = 'K/s'
798                hom(:,2,85,:) = SPREAD( zu, 2, statistic_regions+1 )
799             ENDIF
800
801          CASE ( 'td_nud_q' )
802             IF (  .NOT.  nudging )  THEN
803                message_string = 'data_output_pr = ' //                        &
804                                 TRIM( data_output_pr(var_count) ) //          &
805                                 ' is not implemented for ' //                 &
806                                 'nudging = .FALSE.'
807                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0394',    &
808                               1, 2, 0, 6, 0 )
809             ELSE
810                dopr_index(var_count) = 86
811                dopr_unit             = 'kg/kgs'
812                unit                  = 'kg/kgs'
813                hom(:,2,86,:) = SPREAD( zu, 2, statistic_regions+1 )
814             ENDIF
815
816          CASE ( 'td_nud_u' )
817             IF (  .NOT.  nudging )  THEN
818                message_string = 'data_output_pr = ' //                        &
819                                 TRIM( data_output_pr(var_count) ) //          &
820                                 ' is not implemented for ' //                 &
821                                 'nudging = .FALSE.'
822                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0394',    &
823                               1, 2, 0, 6, 0 )
824             ELSE
825                dopr_index(var_count) = 87
826                dopr_unit             = 'm/s2'
827                unit                  = 'm/s2'
828                hom(:,2,87,:) = SPREAD( zu, 2, statistic_regions+1 )
829             ENDIF
830
831          CASE ( 'td_nud_v' )
832             IF (  .NOT.  nudging )  THEN
833                message_string = 'data_output_pr = ' //                        &
834                                 TRIM( data_output_pr(var_count) ) //          &
835                                 ' is not implemented for ' //                 &
836                                 'nudging = .FALSE.'
837                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0394',    &
838                               1, 2, 0, 6, 0 )
839             ELSE
840                dopr_index(var_count) = 88
841                dopr_unit             = 'm/s2'
842                unit                  = 'm/s2'
843                hom(:,2,88,:) = SPREAD( zu, 2, statistic_regions+1 )
844             ENDIF
845
846
847          CASE DEFAULT
848             unit = 'illegal'
849   
850       END SELECT
851
852    END SUBROUTINE lsf_nudging_check_data_output_pr
853
854!------------------------------------------------------------------------------!
855! Description:
856! ------------
857!> @todo Missing subroutine description.
858!------------------------------------------------------------------------------!
859    SUBROUTINE lsf_nudging_header ( io )
860
861       IMPLICIT NONE
862
863       INTEGER(iwp), INTENT(IN) ::  io !< Unit of the output file
864
865       WRITE ( io, 1 )
866       IF ( large_scale_forcing )  THEN
867          WRITE ( io, 3 )
868          WRITE ( io, 4 )
869
870          IF ( large_scale_subsidence )  THEN
871             IF ( .NOT. use_subsidence_tendencies )  THEN
872                WRITE ( io, 5 )
873             ELSE
874                WRITE ( io, 6 )
875             ENDIF
876          ENDIF
877
878          IF ( bc_pt_b == 'dirichlet' )  THEN
879             WRITE ( io, 12 )
880          ELSEIF ( bc_pt_b == 'neumann' )  THEN
881             WRITE ( io, 13 )
882          ENDIF
883
884          IF ( bc_q_b == 'dirichlet' )  THEN
885             WRITE ( io, 14 )
886          ELSEIF ( bc_q_b == 'neumann' )  THEN
887             WRITE ( io, 15 )
888          ENDIF
889
890          WRITE ( io, 7 )
891          IF ( nudging )  THEN
892             WRITE ( io, 10 )
893          ENDIF
894       ELSE
895          WRITE ( io, 2 )
896          WRITE ( io, 11 )
897       ENDIF
898       IF ( large_scale_subsidence )  THEN
899          WRITE ( io, 8 )
900          WRITE ( io, 9 )
901       ENDIF
902
903
9041 FORMAT (//' Large scale forcing and nudging:'/ &
905              ' -------------------------------'/)
9062 FORMAT (' --> No large scale forcing from external is used (default) ')
9073 FORMAT (' --> Large scale forcing from external file LSF_DATA is used: ')
9084 FORMAT ('     - large scale advection tendencies ')
9095 FORMAT ('     - large scale subsidence velocity w_subs ')
9106 FORMAT ('     - large scale subsidence tendencies ')
9117 FORMAT ('     - and geostrophic wind components ug and vg')
9128 FORMAT (' --> Large-scale vertical motion is used in the ', &
913                  'prognostic equation(s) for')
9149 FORMAT ('     the scalar(s) only')
91510 FORMAT (' --> Nudging is used')
91611 FORMAT (' --> No nudging is used (default) ')
91712 FORMAT ('     - prescribed surface values for temperature')
91813 FORMAT ('     - prescribed surface fluxes for temperature')
91914 FORMAT ('     - prescribed surface values for humidity')
92015 FORMAT ('     - prescribed surface fluxes for humidity')
921
922    END SUBROUTINE lsf_nudging_header
923
924!------------------------------------------------------------------------------!
925! Description:
926! ------------
927!> @todo Missing subroutine description.
928!------------------------------------------------------------------------------!
929    SUBROUTINE lsf_init
930
931       USE control_parameters,                                                 &
932           ONLY:  bc_lr_cyc, bc_ns_cyc
933
934       USE netcdf_data_input_mod,                                              &
935           ONLY:  netcdf_data_input_lsf 
936
937       IMPLICIT NONE
938
939       CHARACTER(100) ::  chmess      !<
940       CHARACTER(1)   ::  hash        !<
941
942       INTEGER(iwp) ::  ierrn         !<
943       INTEGER(iwp) ::  finput = 90   !<
944       INTEGER(iwp) ::  k             !<
945       INTEGER(iwp) ::  nt            !<
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.