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

Last change on this file since 3046 was 3046, checked in by Giersch, 3 years ago

Remaining error messages revised, comments extended

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