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

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

Bugfix in large-scale forcing

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