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

Last change on this file since 2963 was 2938, checked in by suehring, 6 years ago

Nesting in RANS-LES and RANS-RANS mode enabled; synthetic turbulence generator at all lateral boundaries in nesting or non-cyclic forcing mode; revised Inifor initialization in nesting mode

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