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

Last change on this file since 3799 was 3719, checked in by kanani, 5 years ago

Correct and clean-up cpu_logs, some overlapping counts (chemistry_model_mod, disturb_heatflux, large_scale_forcing_nudging_mod, ocean_mod, palm, prognostic_equations, synthetic_turbulence_generator_mod, time_integration, time_integration_spinup, turbulence_closure_mod)

  • Property svn:keywords set to Id
File size: 55.3 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-2019 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: large_scale_forcing_nudging_mod.f90 3719 2019-02-06 13:10:18Z forkel $
27! Removed USE cpulog (unused)
28!
29! 3655 2019-01-07 16:51:22Z knoop
30! unused variables removed
31!
32! 3428 2018-10-25 12:32:05Z gronemeier
33! Rename td_XXX_lpt to td_XXX_thetal
34!
35! 3347 2018-10-15 14:21:08Z suehring
36! Offline nesting parts are moved to an independent module nesting_offl_mod
37!
38! 3341 2018-10-15 10:31:27Z suehring
39! ocean renamed ocean_mode
40!
41! 3241 2018-09-12 15:02:00Z raasch
42! unused variables removed
43!
44! 3183 2018-07-27 14:25:55Z suehring
45! * Adjustment to new Inifor version:
46!   - No vertical interpolation/extrapolation of lateral boundary data required
47!     any more (Inifor can treat grid stretching now
48!   - Revise initialization in case of COSMO forcing
49! * Rename variables and subroutines for offline nesting
50!
51! 3182 2018-07-27 13:36:03Z suehring
52! Error messages revised
53!
54! 3045 2018-05-28 07:55:41Z Giersch
55! Error messages revised
56!
57! 3026 2018-05-22 10:30:53Z schwenkel
58! Changed the name specific humidity to mixing ratio, since we are computing
59! mixing ratios.
60!
61! 2970 2018-04-13 15:09:23Z suehring
62! Bugfix in old large-scale forcing mode
63!
64! 2938 2018-03-27 15:52:42Z suehring
65! Further improvements for nesting in larger-scale model
66!
67! 2863 2018-03-08 11:36:25Z suehring
68! Corrected "Former revisions" section
69!
70! 2696 2017-12-14 17:12:51Z kanani
71! Change in file header (GPL part)
72! Forcing with larger-scale models implemented (MS)
73!
74! 2342 2017-08-08 11:00:43Z boeske
75! fixed check if surface forcing data is available until end of simulation
76!
77! 2320 2017-07-21 12:47:43Z suehring
78! initial revision
79!
80! Description:
81! ------------
82!> Calculates large scale forcings (geostrophic wind and subsidence velocity) as
83!> well as surfaces fluxes dependent on time given in an external file (LSF_DATA).
84!> Moreover, module contains nudging routines, where u, v, pt and q are nudged
85!> to given profiles on a relaxation timescale tnudge.
86!> Profiles are read in from NUDGING_DATA.
87!> Code is based on Neggers et al. (2012) and also in parts on DALES and UCLA-LES.
88!> @todo: Revise reading of ASCII-files
89!> @todo: Remove unused variables and control flags
90!> @todo: Revise large-scale facing of surface variables
91!> @todo: Revise control flags lsf_exception, lsf_surf, lsf_vert, etc.
92!--------------------------------------------------------------------------------!
93 MODULE lsf_nudging_mod
94
95    USE arrays_3d,                                                             &
96        ONLY:  dzw, e, diss, heatflux_input_conversion, pt, pt_init, q,        &
97               q_init, s, tend, u, u_init, ug, v, v_init, vg, w, w_subs,       &
98               waterflux_input_conversion, zu, zw                 
99
100    USE control_parameters,                                                    &
101        ONLY:  bc_lr, bc_ns, bc_pt_b, bc_q_b, constant_diffusion,              &
102               constant_heatflux, constant_waterflux,                          &
103               data_output_pr, dt_3d, end_time,                                &
104               humidity, initializing_actions, intermediate_timestep_count,    &
105               ibc_pt_b, ibc_q_b,                                              &
106               large_scale_forcing, large_scale_subsidence, lsf_surf, lsf_vert,&
107               lsf_exception, message_string, neutral,                         &
108               nudging, passive_scalar, pt_surface, ocean_mode, q_surface,     &
109               surface_heatflux, surface_pressure, surface_waterflux,          &
110               topography, use_subsidence_tendencies
111               
112    USE grid_variables
113
114    USE indices,                                                               &
115        ONLY:  nbgp, ngp_sums_ls, nx, nxl, nxlg, nxlu, nxr, nxrg, ny, nys,     &
116               nysv, nysg, nyn, nyng, nzb, nz, nzt, wall_flags_0
117
118    USE kinds
119
120    USE pegrid
121
122    USE surface_mod,                                                           &
123        ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
124
125    USE statistics,                                                            &
126        ONLY:  hom, statistic_regions, sums_ls_l, weight_substep
127
128    INTEGER(iwp) ::  nlsf = 1000                       !< maximum number of profiles in LSF_DATA (large scale forcing)
129    INTEGER(iwp) ::  ntnudge = 1000                    !< maximum number of profiles in NUDGING_DATA (nudging)
130
131    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ptnudge     !< vertical profile of pot. temperature interpolated to vertical grid (nudging)
132    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  qnudge      !< vertical profile of water vapor mixing ratio interpolated to vertical grid (nudging)
133    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tnudge      !< vertical profile of nudging time scale interpolated to vertical grid (nudging) 
134    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  td_lsa_lpt  !< temperature tendency due to large scale advection (large scale forcing)
135    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  td_lsa_q    !< water vapor mixing ratio tendency due to large scale advection (large scale forcing)
136    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  td_sub_lpt  !< temperature tendency due to subsidence/ascent (large scale forcing)
137    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  td_sub_q    !< water vapor mixing ratio tendency due to subsidence/ascent (large scale forcing)
138    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ug_vert     !< vertical profile of geostrophic wind component in x-direction interpolated to vertical grid (large scale forcing)
139    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  unudge      !< vertical profile of wind component in x-direction interpolated to vertical grid (nudging)
140    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  vnudge      !< vertical profile of wind component in y-direction interpolated to vertical grid (nudging)
141    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  vg_vert     !< vertical profile of geostrophic wind component in y-direction interpolated to vertical grid (large scale forcing)
142    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  wnudge      !< vertical profile of subsidence/ascent velocity interpolated to vertical grid (nudging) ???
143    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  wsubs_vert  !< vertical profile of wind component in z-direction interpolated to vertical grid (nudging) ???
144
145    REAL(wp), DIMENSION(:), ALLOCATABLE ::  shf_surf      !< time-dependent surface sensible heat flux (large scale forcing)
146    REAL(wp), DIMENSION(:), ALLOCATABLE ::  timenudge     !< times at which vertical profiles are defined in NUDGING_DATA (nudging)
147    REAL(wp), DIMENSION(:), ALLOCATABLE ::  time_surf     !< times at which surface values/fluxes are defined in LSF_DATA (large scale forcing)
148    REAL(wp), DIMENSION(:), ALLOCATABLE ::  time_vert     !< times at which vertical profiles are defined in LSF_DATA (large scale forcing)
149
150    REAL(wp), DIMENSION(:), ALLOCATABLE ::  tmp_tnudge    !< current nudging time scale
151
152    REAL(wp), DIMENSION(:), ALLOCATABLE ::  p_surf        !< time-dependent surface pressure (large scale forcing)
153    REAL(wp), DIMENSION(:), ALLOCATABLE ::  pt_surf       !< time-dependent surface temperature (large scale forcing)
154    REAL(wp), DIMENSION(:), ALLOCATABLE ::  qsws_surf     !< time-dependent surface latent heat flux (large scale forcing)
155    REAL(wp), DIMENSION(:), ALLOCATABLE ::  q_surf        !< time-dependent surface water vapor mixing ratio (large scale forcing)
156
157    SAVE
158    PRIVATE
159!
160!-- Public subroutines
161    PUBLIC calc_tnudge, ls_forcing_surf, ls_forcing_vert, ls_advec, lsf_init,  &
162           lsf_nudging_check_parameters, nudge_init,                           &
163           lsf_nudging_check_data_output_pr, lsf_nudging_header,               &
164           nudge, nudge_ref
165           
166!
167!-- Public variables
168    PUBLIC qsws_surf, shf_surf, td_lsa_lpt, td_lsa_q, td_sub_lpt,              &
169           td_sub_q, time_vert
170
171
172    INTERFACE ls_advec
173       MODULE PROCEDURE ls_advec
174       MODULE PROCEDURE ls_advec_ij
175    END INTERFACE ls_advec
176
177    INTERFACE nudge
178       MODULE PROCEDURE nudge
179       MODULE PROCEDURE nudge_ij
180    END INTERFACE nudge
181
182 CONTAINS
183
184
185!------------------------------------------------------------------------------!
186! Description:
187! ------------
188!> @todo Missing subroutine description.
189!------------------------------------------------------------------------------!
190    SUBROUTINE lsf_nudging_check_parameters
191
192       IMPLICIT NONE
193!
194!--    Check nudging and large scale forcing from external file
195       IF ( nudging  .AND.  (  .NOT.  large_scale_forcing ) )  THEN
196          message_string = 'Nudging requires large_scale_forcing = .T.. &'//   &
197                        'Surface fluxes and geostrophic wind should be &'//    &
198                        'prescribed in file LSF_DATA'
199          CALL message( 'check_parameters', 'PA0374', 1, 2, 0, 6, 0 )
200       ENDIF
201
202       IF ( large_scale_forcing  .AND.  ( bc_lr /= 'cyclic'  .OR.              &
203                                          bc_ns /= 'cyclic' ) )  THEN
204          message_string = 'Non-cyclic lateral boundaries do not allow for &'//&
205                        'the usage of large scale forcing from external file.'
206          CALL message( 'check_parameters', 'PA0375', 1, 2, 0, 6, 0 )
207       ENDIF
208
209       IF ( large_scale_forcing  .AND.  (  .NOT.  humidity ) )  THEN
210          message_string = 'The usage of large scale forcing from external &'//&
211                        'file LSF_DATA requires humidity = .T..'
212          CALL message( 'check_parameters', 'PA0376', 1, 2, 0, 6, 0 )
213       ENDIF
214
215       IF ( large_scale_forcing  .AND.  passive_scalar )  THEN
216          message_string = 'The usage of large scale forcing from external &'// &
217                        'file LSF_DATA is not implemented for passive scalars'
218          CALL message( 'check_parameters', 'PA0440', 1, 2, 0, 6, 0 )
219       ENDIF
220
221       IF ( large_scale_forcing  .AND.  topography /= 'flat'                   &
222                              .AND.  .NOT.  lsf_exception )  THEN
223          message_string = 'The usage of large scale forcing from external &'//&
224                        'file LSF_DATA is not implemented for non-flat topography'
225          CALL message( 'check_parameters', 'PA0377', 1, 2, 0, 6, 0 )
226       ENDIF
227
228       IF ( large_scale_forcing  .AND.  ocean_mode )  THEN
229          message_string = 'The usage of large scale forcing from external &'//&
230                        'file LSF_DATA is not implemented for ocean mode'
231          CALL message( 'check_parameters', 'PA0378', 1, 2, 0, 6, 0 )
232       ENDIF
233
234    END SUBROUTINE lsf_nudging_check_parameters
235
236!------------------------------------------------------------------------------!
237! Description:
238! ------------
239!> Check data output of profiles for land surface model
240!------------------------------------------------------------------------------!
241    SUBROUTINE lsf_nudging_check_data_output_pr( variable, var_count, unit,    &
242                                                 dopr_unit )
243 
244       USE profil_parameter
245
246       IMPLICIT NONE
247   
248       CHARACTER (LEN=*) ::  unit      !<
249       CHARACTER (LEN=*) ::  variable  !<
250       CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
251 
252       INTEGER(iwp) ::  var_count     !<
253
254       SELECT CASE ( TRIM( variable ) )
255       
256
257          CASE ( 'td_lsa_thetal' )
258             IF (  .NOT.  large_scale_forcing )  THEN
259                message_string = 'data_output_pr = ' //                        &
260                                 TRIM( data_output_pr(var_count) ) //          &
261                                 ' is not implemented for ' //                 &
262                                 'large_scale_forcing = .FALSE.'
263                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393',    &
264                               1, 2, 0, 6, 0 )
265             ELSE
266                dopr_index(var_count) = 81
267                dopr_unit             = 'K/s'
268                unit                  = 'K/s'
269                hom(:,2,81,:) = SPREAD( zu, 2, statistic_regions+1 )
270             ENDIF
271
272          CASE ( 'td_lsa_q' )
273             IF (  .NOT.  large_scale_forcing )  THEN
274                message_string = 'data_output_pr = ' //                        &
275                                 TRIM( data_output_pr(var_count) ) //          &
276                                 ' is not implemented for ' //                 &
277                                 'large_scale_forcing = .FALSE.'
278                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393',    &
279                               1, 2, 0, 6, 0 )
280             ELSE
281                dopr_index(var_count) = 82
282                dopr_unit             = 'kg/kgs'
283                unit                  = 'kg/kgs'
284                hom(:,2,82,:) = SPREAD( zu, 2, statistic_regions+1 )
285             ENDIF
286          CASE ( 'td_sub_thetal' )
287             IF (  .NOT.  large_scale_forcing )  THEN
288                message_string = 'data_output_pr = ' //                        &
289                                 TRIM( data_output_pr(var_count) ) //          &
290                                 ' is not implemented for ' //                 &
291                                 'large_scale_forcing = .FALSE.'
292                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393',    &
293                               1, 2, 0, 6, 0 )
294             ELSE
295                dopr_index(var_count) = 83
296                dopr_unit             = 'K/s'
297                unit                  = 'K/s'
298                hom(:,2,83,:) = SPREAD( zu, 2, statistic_regions+1 )
299             ENDIF
300
301          CASE ( 'td_sub_q' )
302             IF (  .NOT.  large_scale_forcing )  THEN
303                message_string = 'data_output_pr = ' //                        &
304                                 TRIM( data_output_pr(var_count) ) //          &
305                                 ' is not implemented for ' //                 &
306                                 'large_scale_forcing = .FALSE.'
307                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393',    &
308                               1, 2, 0, 6, 0 )
309             ELSE
310                dopr_index(var_count) = 84
311                dopr_unit             = 'kg/kgs'
312                unit                  = 'kg/kgs'
313                hom(:,2,84,:) = SPREAD( zu, 2, statistic_regions+1 )
314             ENDIF
315
316          CASE ( 'td_nud_thetal' )
317             IF (  .NOT.  nudging )  THEN
318                message_string = 'data_output_pr = ' //                        &
319                                 TRIM( data_output_pr(var_count) ) //          &
320                                 ' is not implemented for ' //                 &
321                                 'nudging = .FALSE.'
322                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0394',    &
323                               1, 2, 0, 6, 0 )
324             ELSE
325                dopr_index(var_count) = 85
326                dopr_unit             = 'K/s'
327                unit                  = 'K/s'
328                hom(:,2,85,:) = SPREAD( zu, 2, statistic_regions+1 )
329             ENDIF
330
331          CASE ( 'td_nud_q' )
332             IF (  .NOT.  nudging )  THEN
333                message_string = 'data_output_pr = ' //                        &
334                                 TRIM( data_output_pr(var_count) ) //          &
335                                 ' is not implemented for ' //                 &
336                                 'nudging = .FALSE.'
337                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0394',    &
338                               1, 2, 0, 6, 0 )
339             ELSE
340                dopr_index(var_count) = 86
341                dopr_unit             = 'kg/kgs'
342                unit                  = 'kg/kgs'
343                hom(:,2,86,:) = SPREAD( zu, 2, statistic_regions+1 )
344             ENDIF
345
346          CASE ( 'td_nud_u' )
347             IF (  .NOT.  nudging )  THEN
348                message_string = 'data_output_pr = ' //                        &
349                                 TRIM( data_output_pr(var_count) ) //          &
350                                 ' is not implemented for ' //                 &
351                                 'nudging = .FALSE.'
352                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0394',    &
353                               1, 2, 0, 6, 0 )
354             ELSE
355                dopr_index(var_count) = 87
356                dopr_unit             = 'm/s2'
357                unit                  = 'm/s2'
358                hom(:,2,87,:) = SPREAD( zu, 2, statistic_regions+1 )
359             ENDIF
360
361          CASE ( 'td_nud_v' )
362             IF (  .NOT.  nudging )  THEN
363                message_string = 'data_output_pr = ' //                        &
364                                 TRIM( data_output_pr(var_count) ) //          &
365                                 ' is not implemented for ' //                 &
366                                 'nudging = .FALSE.'
367                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0394',    &
368                               1, 2, 0, 6, 0 )
369             ELSE
370                dopr_index(var_count) = 88
371                dopr_unit             = 'm/s2'
372                unit                  = 'm/s2'
373                hom(:,2,88,:) = SPREAD( zu, 2, statistic_regions+1 )
374             ENDIF
375
376
377          CASE DEFAULT
378             unit = 'illegal'
379   
380       END SELECT
381
382    END SUBROUTINE lsf_nudging_check_data_output_pr
383
384!------------------------------------------------------------------------------!
385! Description:
386! ------------
387!> @todo Missing subroutine description.
388!------------------------------------------------------------------------------!
389    SUBROUTINE lsf_nudging_header ( io )
390
391       IMPLICIT NONE
392
393       INTEGER(iwp), INTENT(IN) ::  io !< Unit of the output file
394
395       WRITE ( io, 1 )
396       IF ( large_scale_forcing )  THEN
397          WRITE ( io, 3 )
398          WRITE ( io, 4 )
399
400          IF ( large_scale_subsidence )  THEN
401             IF ( .NOT. use_subsidence_tendencies )  THEN
402                WRITE ( io, 5 )
403             ELSE
404                WRITE ( io, 6 )
405             ENDIF
406          ENDIF
407
408          IF ( bc_pt_b == 'dirichlet' )  THEN
409             WRITE ( io, 12 )
410          ELSEIF ( bc_pt_b == 'neumann' )  THEN
411             WRITE ( io, 13 )
412          ENDIF
413
414          IF ( bc_q_b == 'dirichlet' )  THEN
415             WRITE ( io, 14 )
416          ELSEIF ( bc_q_b == 'neumann' )  THEN
417             WRITE ( io, 15 )
418          ENDIF
419
420          WRITE ( io, 7 )
421          IF ( nudging )  THEN
422             WRITE ( io, 10 )
423          ENDIF
424       ELSE
425          WRITE ( io, 2 )
426          WRITE ( io, 11 )
427       ENDIF
428       IF ( large_scale_subsidence )  THEN
429          WRITE ( io, 8 )
430          WRITE ( io, 9 )
431       ENDIF
432
433
4341 FORMAT (//' Large scale forcing and nudging:'/ &
435              ' -------------------------------'/)
4362 FORMAT (' --> No large scale forcing from external is used (default) ')
4373 FORMAT (' --> Large scale forcing from external file LSF_DATA is used: ')
4384 FORMAT ('     - large scale advection tendencies ')
4395 FORMAT ('     - large scale subsidence velocity w_subs ')
4406 FORMAT ('     - large scale subsidence tendencies ')
4417 FORMAT ('     - and geostrophic wind components ug and vg')
4428 FORMAT (' --> Large-scale vertical motion is used in the ', &
443                  'prognostic equation(s) for')
4449 FORMAT ('     the scalar(s) only')
44510 FORMAT (' --> Nudging is used')
44611 FORMAT (' --> No nudging is used (default) ')
44712 FORMAT ('     - prescribed surface values for temperature')
44813 FORMAT ('     - prescribed surface fluxes for temperature')
44914 FORMAT ('     - prescribed surface values for humidity')
45015 FORMAT ('     - prescribed surface fluxes for humidity')
451
452    END SUBROUTINE lsf_nudging_header 
453
454!------------------------------------------------------------------------------!
455! Description:
456! ------------
457!> @todo Missing subroutine description.
458!------------------------------------------------------------------------------!
459    SUBROUTINE lsf_init
460
461       IMPLICIT NONE
462
463       CHARACTER(100) ::  chmess      !<
464       CHARACTER(1)   ::  hash        !<
465
466       INTEGER(iwp) ::  ierrn         !<
467       INTEGER(iwp) ::  finput = 90   !<
468       INTEGER(iwp) ::  k             !<
469       INTEGER(iwp) ::  nt            !<
470
471       REAL(wp) ::  fac               !<
472       REAL(wp) ::  highheight        !<
473       REAL(wp) ::  highug_vert       !<
474       REAL(wp) ::  highvg_vert       !<
475       REAL(wp) ::  highwsubs_vert    !<
476       REAL(wp) ::  lowheight         !<
477       REAL(wp) ::  lowug_vert        !<
478       REAL(wp) ::  lowvg_vert        !<
479       REAL(wp) ::  lowwsubs_vert     !<
480       REAL(wp) ::  high_td_lsa_lpt   !<
481       REAL(wp) ::  low_td_lsa_lpt    !<
482       REAL(wp) ::  high_td_lsa_q     !<
483       REAL(wp) ::  low_td_lsa_q      !<
484       REAL(wp) ::  high_td_sub_lpt   !<
485       REAL(wp) ::  low_td_sub_lpt    !<
486       REAL(wp) ::  high_td_sub_q     !<
487       REAL(wp) ::  low_td_sub_q      !<
488       REAL(wp) ::  r_dummy           !<
489
490       ALLOCATE( p_surf(0:nlsf), pt_surf(0:nlsf), q_surf(0:nlsf),              &
491                 qsws_surf(0:nlsf), shf_surf(0:nlsf),                          &
492                 td_lsa_lpt(nzb:nzt+1,0:nlsf), td_lsa_q(nzb:nzt+1,0:nlsf),     &
493                 td_sub_lpt(nzb:nzt+1,0:nlsf), td_sub_q(nzb:nzt+1,0:nlsf),     &
494                 time_vert(0:nlsf), time_surf(0:nlsf),                         &
495                 ug_vert(nzb:nzt+1,0:nlsf), vg_vert(nzb:nzt+1,0:nlsf),         &
496                 wsubs_vert(nzb:nzt+1,0:nlsf) )
497
498       p_surf = 0.0_wp; pt_surf = 0.0_wp; q_surf = 0.0_wp; qsws_surf = 0.0_wp
499       shf_surf = 0.0_wp; time_vert = 0.0_wp; td_lsa_lpt = 0.0_wp
500       td_lsa_q = 0.0_wp; td_sub_lpt = 0.0_wp; td_sub_q = 0.0_wp
501       time_surf = 0.0_wp; ug_vert = 0.0_wp; vg_vert = 0.0_wp
502       wsubs_vert = 0.0_wp
503
504!
505!--    Array for storing large scale forcing and nudging tendencies at each
506!--    timestep for data output
507       ALLOCATE( sums_ls_l(nzb:nzt+1,0:7) )
508       sums_ls_l = 0.0_wp
509
510       ngp_sums_ls = (nz+2)*6
511
512       OPEN ( finput, FILE='LSF_DATA', STATUS='OLD', &
513              FORM='FORMATTED', IOSTAT=ierrn )
514
515       IF ( ierrn /= 0 )  THEN
516          message_string = 'file LSF_DATA does not exist'
517          CALL message( 'ls_forcing', 'PA0368', 1, 2, 0, 6, 0 )
518       ENDIF
519
520       ierrn = 0
521!
522!--    First three lines of LSF_DATA contain header
523       READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
524       READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
525       READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
526
527       IF ( ierrn /= 0 )  THEN
528          message_string = 'errors in file LSF_DATA'
529          CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
530       ENDIF
531
532!
533!--    Surface values are read in
534       nt     = 0
535       ierrn = 0
536
537       DO WHILE ( time_surf(nt) < end_time )
538          nt = nt + 1
539          READ ( finput, *, IOSTAT = ierrn ) time_surf(nt), shf_surf(nt),      &
540                                             qsws_surf(nt), pt_surf(nt),       &
541                                             q_surf(nt), p_surf(nt)
542
543          IF ( ierrn /= 0 )  THEN
544            WRITE ( message_string, * ) 'No time dependent surface ' //        &
545                              'variables in & LSF_DATA for end of run found'
546
547             CALL message( 'ls_forcing', 'PA0363', 1, 2, 0, 6, 0 )
548          ENDIF
549       ENDDO
550
551       IF ( time_surf(1) > end_time )  THEN
552          WRITE ( message_string, * ) 'Time dependent surface variables in ' //&
553                                      '&LSF_DATA set in after end of ' ,       &
554                                      'simulation - lsf_surf is set to FALSE'
555          CALL message( 'ls_forcing', 'PA0371', 0, 0, 0, 6, 0 )
556          lsf_surf = .FALSE.
557       ENDIF
558
559!
560!--    Go to the end of the list with surface variables
561       DO WHILE ( ierrn == 0 )
562          READ ( finput, *, IOSTAT = ierrn ) r_dummy
563       ENDDO
564
565!
566!--    Profiles of ug, vg and w_subs are read in (large scale forcing)
567
568       nt = 0
569       DO WHILE ( time_vert(nt) < end_time )
570          nt = nt + 1
571          hash = "#"
572          ierrn = 1 ! not zero
573!
574!--       Search for the next line consisting of "# time",
575!--       from there onwards the profiles will be read
576          DO WHILE ( .NOT. ( hash == "#" .AND. ierrn == 0 ) ) 
577             READ ( finput, *, IOSTAT=ierrn ) hash, time_vert(nt)
578             IF ( ierrn < 0 )  THEN
579                WRITE( message_string, * ) 'No time dependent vertical profiles',&
580                                 ' in & LSF_DATA for end of run found'
581                CALL message( 'ls_forcing', 'PA0372', 1, 2, 0, 6, 0 )
582             ENDIF
583          ENDDO
584
585          IF ( nt == 1 .AND. time_vert(nt) > end_time ) EXIT
586
587          READ ( finput, *, IOSTAT=ierrn ) lowheight, lowug_vert, lowvg_vert,  &
588                                           lowwsubs_vert, low_td_lsa_lpt,      &
589                                           low_td_lsa_q, low_td_sub_lpt,       &
590                                           low_td_sub_q
591          IF ( ierrn /= 0 )  THEN
592             message_string = 'errors in file LSF_DATA'
593             CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
594          ENDIF
595
596          READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert,            &
597                                           highvg_vert, highwsubs_vert,        &
598                                           high_td_lsa_lpt, high_td_lsa_q,     &
599                                           high_td_sub_lpt, high_td_sub_q
600       
601          IF ( ierrn /= 0 )  THEN
602             message_string = 'errors in file LSF_DATA'
603             CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
604          ENDIF
605
606
607          DO  k = nzb, nzt+1
608             IF ( highheight < zu(k) )  THEN
609                lowheight      = highheight
610                lowug_vert     = highug_vert
611                lowvg_vert     = highvg_vert
612                lowwsubs_vert  = highwsubs_vert
613                low_td_lsa_lpt = high_td_lsa_lpt
614                low_td_lsa_q   = high_td_lsa_q
615                low_td_sub_lpt = high_td_sub_lpt
616                low_td_sub_q   = high_td_sub_q
617
618                ierrn = 0
619                READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert,      &
620                                                 highvg_vert, highwsubs_vert,  &
621                                                 high_td_lsa_lpt,              &
622                                                 high_td_lsa_q,                &
623                                                 high_td_sub_lpt, high_td_sub_q
624
625                IF ( ierrn /= 0 )  THEN
626                   WRITE( message_string, * ) 'zu(',k,') = ', zu(k), 'm ',     &
627                        'is higher than the maximum height in LSF_DATA ',      &
628                        'which is ', lowheight, 'm. Interpolation on PALM ',   &
629                        'grid is not possible.'
630                   CALL message( 'ls_forcing', 'PA0395', 1, 2, 0, 6, 0 )
631                ENDIF
632
633             ENDIF
634
635!
636!--          Interpolation of prescribed profiles in space
637             fac = (highheight-zu(k))/(highheight - lowheight)
638
639             ug_vert(k,nt)    = fac * lowug_vert                               &
640                                + ( 1.0_wp - fac ) * highug_vert
641             vg_vert(k,nt)    = fac * lowvg_vert                               &
642                                + ( 1.0_wp - fac ) * highvg_vert
643             wsubs_vert(k,nt) = fac * lowwsubs_vert                            &
644                                + ( 1.0_wp - fac ) * highwsubs_vert
645
646             td_lsa_lpt(k,nt) = fac * low_td_lsa_lpt                           &
647                                + ( 1.0_wp - fac ) * high_td_lsa_lpt
648             td_lsa_q(k,nt)   = fac * low_td_lsa_q                             &
649                                + ( 1.0_wp - fac ) * high_td_lsa_q
650             td_sub_lpt(k,nt) = fac * low_td_sub_lpt                           &
651                                + ( 1.0_wp - fac ) * high_td_sub_lpt
652             td_sub_q(k,nt)   = fac * low_td_sub_q                             &
653                                + ( 1.0_wp - fac ) * high_td_sub_q
654
655          ENDDO
656
657       ENDDO 
658
659!
660!--    Large scale vertical velocity has to be zero at the surface
661       wsubs_vert(nzb,:) = 0.0_wp
662   
663       IF ( time_vert(1) > end_time )  THEN
664          WRITE ( message_string, * ) 'Time dependent large scale profile ',   &
665                             'forcing from&LSF_DATA sets in after end of ' ,   &
666                             'simulation - lsf_vert is set to FALSE'
667          CALL message( 'ls_forcing', 'PA0373', 0, 0, 0, 6, 0 )
668          lsf_vert = .FALSE.
669       ENDIF
670
671       CLOSE( finput )
672
673    END SUBROUTINE lsf_init
674
675!------------------------------------------------------------------------------!
676! Description:
677! ------------
678!> @todo Missing subroutine description.
679!------------------------------------------------------------------------------!
680    SUBROUTINE ls_forcing_surf ( time )
681
682       IMPLICIT NONE
683
684       INTEGER(iwp) ::  nt                     !<
685
686       REAL(wp)             :: dum_surf_flux  !<
687       REAL(wp)             :: fac            !<
688       REAL(wp), INTENT(in) :: time           !<
689
690!
691!--    Interpolation in time of LSF_DATA at the surface
692       nt = 1
693       DO WHILE ( time > time_surf(nt) )
694          nt = nt + 1
695       ENDDO
696       IF ( time /= time_surf(nt) )  THEN
697         nt = nt - 1
698       ENDIF
699
700       fac = ( time -time_surf(nt) ) / ( time_surf(nt+1) - time_surf(nt) )
701
702       IF ( ibc_pt_b == 0 )  THEN
703!
704!--       In case of Dirichlet boundary condition shf must not
705!--       be set - it is calculated via MOST in prandtl_fluxes
706          pt_surface = pt_surf(nt) + fac * ( pt_surf(nt+1) - pt_surf(nt) )
707
708       ELSEIF ( ibc_pt_b == 1 )  THEN
709!
710!--       In case of Neumann boundary condition pt_surface is needed for
711!--       calculation of reference density
712          dum_surf_flux = ( shf_surf(nt) + fac *                               &
713                            ( shf_surf(nt+1) - shf_surf(nt) )                  &
714                          ) * heatflux_input_conversion(nzb)
715!
716!--       Save surface sensible heat flux on default, natural and urban surface
717!--       type, if required
718          IF ( surf_def_h(0)%ns >= 1 )  surf_def_h(0)%shf(:) = dum_surf_flux
719          IF ( surf_lsm_h%ns    >= 1 )  surf_lsm_h%shf(:)    = dum_surf_flux
720          IF ( surf_usm_h%ns    >= 1 )  surf_usm_h%shf(:)    = dum_surf_flux
721
722          pt_surface    = pt_surf(nt) + fac * ( pt_surf(nt+1) - pt_surf(nt) )
723
724       ENDIF
725
726       IF ( ibc_q_b == 0 )  THEN
727!
728!--       In case of Dirichlet boundary condition qsws must not
729!--       be set - it is calculated via MOST in prandtl_fluxes
730          q_surface = q_surf(nt) + fac * ( q_surf(nt+1) - q_surf(nt) )
731
732       ELSEIF ( ibc_q_b == 1 )  THEN
733          dum_surf_flux = ( qsws_surf(nt) + fac *                              &
734                             ( qsws_surf(nt+1) - qsws_surf(nt) )               &
735                             ) * waterflux_input_conversion(nzb)
736!
737!--       Save surface sensible heat flux on default, natural and urban surface
738!--       type, if required
739          IF ( surf_def_h(0)%ns >= 1 )  surf_def_h(0)%qsws(:) = dum_surf_flux
740          IF ( surf_lsm_h%ns    >= 1 )  surf_lsm_h%qsws(:)    = dum_surf_flux
741          IF ( surf_usm_h%ns    >= 1 )  surf_usm_h%qsws(:)    = dum_surf_flux
742
743       ENDIF
744!
745!--    Surface heat- and waterflux will be written later onto surface elements
746       IF ( .NOT.  neutral  .AND.  constant_heatflux  .AND.                    &
747            TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
748             surface_heatflux = shf_surf(1)
749       ENDIF
750       IF ( humidity  .AND.  constant_waterflux  .AND.                         &
751            TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
752             surface_waterflux = qsws_surf(1)
753       ENDIF
754
755       surface_pressure = p_surf(nt) + fac * ( p_surf(nt+1) - p_surf(nt) )
756
757    END SUBROUTINE ls_forcing_surf 
758
759
760
761
762!------------------------------------------------------------------------------!
763! Description:
764! ------------
765!> @todo Missing subroutine description.
766!------------------------------------------------------------------------------!
767    SUBROUTINE ls_forcing_vert ( time )
768
769
770       IMPLICIT NONE
771
772       INTEGER(iwp) ::  nt                     !<
773
774       REAL(wp)             ::  fac           !<
775       REAL(wp), INTENT(in) ::  time          !<
776
777!
778!--    Interpolation in time of LSF_DATA for ug, vg and w_subs
779       nt = 1
780       DO WHILE ( time > time_vert(nt) )
781          nt = nt + 1
782       ENDDO
783       IF ( time /= time_vert(nt) )  THEN
784         nt = nt - 1
785       ENDIF
786
787       fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
788
789       ug     = ug_vert(:,nt) + fac * ( ug_vert(:,nt+1) - ug_vert(:,nt) )
790       vg     = vg_vert(:,nt) + fac * ( vg_vert(:,nt+1) - vg_vert(:,nt) )
791
792       IF ( large_scale_subsidence )  THEN
793          w_subs = wsubs_vert(:,nt)                                            &
794                   + fac * ( wsubs_vert(:,nt+1) - wsubs_vert(:,nt) )
795       ENDIF
796
797    END SUBROUTINE ls_forcing_vert
798
799
800!------------------------------------------------------------------------------!
801! Description:
802! ------------
803!> Call for all grid points
804!------------------------------------------------------------------------------!
805    SUBROUTINE ls_advec ( time, prog_var )
806     
807
808       IMPLICIT NONE
809
810       CHARACTER (LEN=*) ::  prog_var   !<
811
812       REAL(wp), INTENT(in)  :: time    !<
813       REAL(wp) :: fac                  !< 
814
815       INTEGER(iwp) ::  i               !<
816       INTEGER(iwp) ::  j               !<
817       INTEGER(iwp) ::  k               !<
818       INTEGER(iwp) ::  nt               !<
819
820!
821!--    Interpolation in time of LSF_DATA
822       nt = 1
823       DO WHILE ( time > time_vert(nt) )
824          nt = nt + 1
825       ENDDO
826       IF ( time /= time_vert(nt) )  THEN
827         nt = nt - 1
828       ENDIF
829
830       fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
831
832!
833!--    Add horizontal large scale advection tendencies of pt and q
834       SELECT CASE ( prog_var )
835
836          CASE ( 'pt' )
837
838             DO  i = nxl, nxr
839                DO  j = nys, nyn
840                   DO  k = nzb+1, nzt
841                      tend(k,j,i) = tend(k,j,i) + td_lsa_lpt(k,nt) + fac *     &
842                                    ( td_lsa_lpt(k,nt+1) - td_lsa_lpt(k,nt) ) *&
843                                        MERGE( 1.0_wp, 0.0_wp,                 &
844                                               BTEST( wall_flags_0(k,j,i), 0 ) )
845                   ENDDO
846                ENDDO
847             ENDDO
848
849          CASE ( 'q' )
850
851             DO  i = nxl, nxr
852                DO  j = nys, nyn
853                   DO  k = nzb+1, nzt
854                      tend(k,j,i) = tend(k,j,i) + td_lsa_q(k,nt) + fac *       &
855                                    ( td_lsa_q(k,nt+1) - td_lsa_q(k,nt) ) *    &
856                                        MERGE( 1.0_wp, 0.0_wp,                 &
857                                               BTEST( wall_flags_0(k,j,i), 0 ) )
858                   ENDDO
859                ENDDO
860             ENDDO
861
862       END SELECT
863
864!
865!--    Subsidence of pt and q with prescribed subsidence tendencies
866       IF ( large_scale_subsidence .AND. use_subsidence_tendencies )  THEN
867
868          SELECT CASE ( prog_var )
869
870             CASE ( 'pt' )
871
872                DO  i = nxl, nxr
873                   DO  j = nys, nyn
874                      DO  k = nzb+1, nzt
875                         tend(k,j,i) = tend(k,j,i) + td_sub_lpt(k,nt) + fac *  &
876                                     ( td_sub_lpt(k,nt+1) - td_sub_lpt(k,nt) )*&
877                                        MERGE( 1.0_wp, 0.0_wp,                 &
878                                               BTEST( wall_flags_0(k,j,i), 0 ) )
879                      ENDDO
880                   ENDDO
881                ENDDO
882 
883             CASE ( 'q' )
884
885                DO  i = nxl, nxr
886                   DO  j = nys, nyn
887                      DO  k = nzb+1, nzt
888                         tend(k,j,i) = tend(k,j,i) + td_sub_q(k,nt) + fac *    &
889                                       ( td_sub_q(k,nt+1) - td_sub_q(k,nt) ) * &
890                                        MERGE( 1.0_wp, 0.0_wp,                 &
891                                               BTEST( wall_flags_0(k,j,i), 0 ) )
892                      ENDDO
893                   ENDDO
894                ENDDO
895
896          END SELECT
897
898       ENDIF
899
900    END SUBROUTINE ls_advec
901
902
903!------------------------------------------------------------------------------!
904! Description:
905! ------------
906!> Call for grid point i,j
907!------------------------------------------------------------------------------!
908    SUBROUTINE ls_advec_ij ( i, j, time, prog_var )
909
910       IMPLICIT NONE
911
912       CHARACTER (LEN=*) ::  prog_var   !<
913
914       REAL(wp), INTENT(in)  :: time    !<
915       REAL(wp) :: fac                  !<
916
917       INTEGER(iwp) ::  i               !<
918       INTEGER(iwp) ::  j               !<
919       INTEGER(iwp) ::  k               !<
920       INTEGER(iwp) ::  nt               !<
921
922!
923!--    Interpolation in time of LSF_DATA
924       nt = 1
925       DO WHILE ( time > time_vert(nt) )
926          nt = nt + 1
927       ENDDO
928       IF ( time /= time_vert(nt) )  THEN
929         nt = nt - 1
930       ENDIF
931
932       fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
933
934!
935!--    Add horizontal large scale advection tendencies of pt and q
936       SELECT CASE ( prog_var )
937
938          CASE ( 'pt' )
939
940             DO  k = nzb+1, nzt
941                tend(k,j,i) = tend(k,j,i) + td_lsa_lpt(k,nt)                   &
942                             + fac * ( td_lsa_lpt(k,nt+1) - td_lsa_lpt(k,nt) )*&
943                                        MERGE( 1.0_wp, 0.0_wp,                 &
944                                               BTEST( wall_flags_0(k,j,i), 0 ) )
945             ENDDO
946
947          CASE ( 'q' )
948
949             DO  k = nzb+1, nzt
950                tend(k,j,i) = tend(k,j,i) + td_lsa_q(k,nt)                     &
951                              + fac * ( td_lsa_q(k,nt+1) - td_lsa_q(k,nt) ) *  &
952                                        MERGE( 1.0_wp, 0.0_wp,                 &
953                                               BTEST( wall_flags_0(k,j,i), 0 ) )
954             ENDDO
955
956       END SELECT
957
958!
959!--    Subsidence of pt and q with prescribed profiles
960       IF ( large_scale_subsidence .AND. use_subsidence_tendencies )  THEN
961
962          SELECT CASE ( prog_var )
963
964             CASE ( 'pt' )
965
966                DO  k = nzb+1, nzt
967                   tend(k,j,i) = tend(k,j,i) + td_sub_lpt(k,nt) + fac *        &
968                                 ( td_sub_lpt(k,nt+1) - td_sub_lpt(k,nt) ) *   &
969                                        MERGE( 1.0_wp, 0.0_wp,                 &
970                                               BTEST( wall_flags_0(k,j,i), 0 ) )
971                ENDDO
972 
973             CASE ( 'q' )
974
975                DO  k = nzb+1, nzt
976                   tend(k,j,i) = tend(k,j,i) + td_sub_q(k,nt) + fac *          &
977                                 ( td_sub_q(k,nt+1) - td_sub_q(k,nt) ) *       &
978                                        MERGE( 1.0_wp, 0.0_wp,                 &
979                                               BTEST( wall_flags_0(k,j,i), 0 ) )
980                ENDDO
981
982          END SELECT
983
984       ENDIF
985
986    END SUBROUTINE ls_advec_ij
987
988
989!------------------------------------------------------------------------------!
990! Description:
991! ------------
992!> @todo Missing subroutine description.
993!------------------------------------------------------------------------------!
994    SUBROUTINE nudge_init
995
996       IMPLICIT NONE
997
998
999       INTEGER(iwp) ::  finput = 90  !<
1000       INTEGER(iwp) ::  ierrn        !<
1001       INTEGER(iwp) ::  k            !<
1002       INTEGER(iwp) ::  nt            !<
1003
1004       CHARACTER(1) ::  hash     !<
1005
1006       REAL(wp) ::  highheight   !<
1007       REAL(wp) ::  highqnudge   !<
1008       REAL(wp) ::  highptnudge  !<
1009       REAL(wp) ::  highunudge   !<
1010       REAL(wp) ::  highvnudge   !<
1011       REAL(wp) ::  highwnudge   !<
1012       REAL(wp) ::  hightnudge   !<
1013
1014       REAL(wp) ::  lowheight    !<
1015       REAL(wp) ::  lowqnudge    !<
1016       REAL(wp) ::  lowptnudge   !<
1017       REAL(wp) ::  lowunudge    !<
1018       REAL(wp) ::  lowvnudge    !<
1019       REAL(wp) ::  lowwnudge    !<
1020       REAL(wp) ::  lowtnudge    !<
1021
1022       REAL(wp) ::  fac          !<
1023
1024       ALLOCATE( ptnudge(nzb:nzt+1,1:ntnudge), qnudge(nzb:nzt+1,1:ntnudge), &
1025                 tnudge(nzb:nzt+1,1:ntnudge), unudge(nzb:nzt+1,1:ntnudge),  &
1026                 vnudge(nzb:nzt+1,1:ntnudge), wnudge(nzb:nzt+1,1:ntnudge)  )
1027
1028       ALLOCATE( tmp_tnudge(nzb:nzt) )
1029
1030       ALLOCATE( timenudge(0:ntnudge) )
1031
1032       ptnudge = 0.0_wp; qnudge = 0.0_wp; tnudge = 0.0_wp; unudge = 0.0_wp
1033       vnudge = 0.0_wp; wnudge = 0.0_wp; timenudge = 0.0_wp
1034!
1035!--    Initialize array tmp_nudge with a current nudging time scale of 6 hours
1036       tmp_tnudge = 21600.0_wp
1037
1038       nt = 0
1039       OPEN ( finput, FILE='NUDGING_DATA', STATUS='OLD', &
1040              FORM='FORMATTED', IOSTAT=ierrn )
1041
1042       IF ( ierrn /= 0 )  THEN
1043          message_string = 'file NUDGING_DATA does not exist'
1044          CALL message( 'nudging', 'PA0365', 1, 2, 0, 6, 0 )
1045       ENDIF
1046
1047       ierrn = 0
1048
1049 rloop:DO
1050          nt = nt + 1
1051          hash = "#"
1052          ierrn = 1 ! not zero
1053!
1054!--       Search for the next line consisting of "# time",
1055!--       from there onwards the profiles will be read
1056          DO WHILE ( .NOT. ( hash == "#" .AND. ierrn == 0 ) ) 
1057         
1058            READ ( finput, *, IOSTAT=ierrn ) hash, timenudge(nt)
1059            IF ( ierrn < 0 )  EXIT rloop
1060
1061          ENDDO
1062
1063          ierrn = 0
1064          READ ( finput, *, IOSTAT=ierrn ) lowheight, lowtnudge, lowunudge,   &
1065                                           lowvnudge, lowwnudge , lowptnudge, &
1066                                           lowqnudge
1067
1068          IF ( ierrn /= 0 )  THEN
1069             message_string = 'errors in file NUDGING_DATA'
1070             CALL message( 'nudging', 'PA0366', 1, 2, 0, 6, 0 )
1071          ENDIF
1072
1073          ierrn = 0
1074          READ ( finput, *, IOSTAT=ierrn ) highheight, hightnudge, highunudge,   &
1075                                           highvnudge, highwnudge , highptnudge, &
1076                                           highqnudge
1077
1078          IF ( ierrn /= 0 )  THEN
1079             message_string = 'errors in file NUDGING_DATA'
1080             CALL message( 'nudging', 'PA0366', 1, 2, 0, 6, 0 )
1081          ENDIF
1082
1083          DO  k = nzb, nzt+1
1084             DO WHILE ( highheight < zu(k) )
1085                lowheight  = highheight
1086                lowtnudge  = hightnudge
1087                lowunudge  = highunudge
1088                lowvnudge  = highvnudge
1089                lowwnudge  = highwnudge
1090                lowptnudge = highptnudge
1091                lowqnudge  = highqnudge
1092 
1093                ierrn = 0
1094                READ ( finput, *, IOSTAT=ierrn )  highheight , hightnudge ,    &
1095                                                  highunudge , highvnudge ,    &
1096                                                  highwnudge , highptnudge,    &
1097                                                  highqnudge
1098                IF (ierrn /= 0 )  THEN
1099                   WRITE( message_string, * ) 'zu(',k,') = ', zu(k), 'm is ',  &
1100                        'higher than the maximum height in NUDING_DATA which ',&
1101                        'is ', lowheight, 'm. Interpolation on PALM ',         &
1102                        'grid is not possible.'
1103                   CALL message( 'nudging', 'PA0364', 1, 2, 0, 6, 0 )
1104                ENDIF
1105             ENDDO
1106
1107!
1108!--          Interpolation of prescribed profiles in space
1109
1110             fac = ( highheight - zu(k) ) / ( highheight - lowheight )
1111
1112             tnudge(k,nt)  = fac * lowtnudge  + ( 1.0_wp - fac ) * hightnudge
1113             unudge(k,nt)  = fac * lowunudge  + ( 1.0_wp - fac ) * highunudge
1114             vnudge(k,nt)  = fac * lowvnudge  + ( 1.0_wp - fac ) * highvnudge
1115             wnudge(k,nt)  = fac * lowwnudge  + ( 1.0_wp - fac ) * highwnudge
1116             ptnudge(k,nt) = fac * lowptnudge + ( 1.0_wp - fac ) * highptnudge
1117             qnudge(k,nt)  = fac * lowqnudge  + ( 1.0_wp - fac ) * highqnudge
1118          ENDDO
1119
1120       ENDDO rloop
1121
1122       CLOSE ( finput )
1123
1124!
1125!--    Overwrite initial profiles in case of nudging
1126       IF ( nudging )  THEN
1127          pt_init = ptnudge(:,1)
1128          u_init  = unudge(:,1)
1129          v_init  = vnudge(:,1)
1130          IF ( humidity  )  THEN ! is passive_scalar correct???
1131             q_init = qnudge(:,1)
1132          ENDIF
1133
1134          WRITE( message_string, * ) 'Initial profiles of u, v, pt and q ',    &
1135                                     'from NUDGING_DATA are used.'
1136          CALL message( 'large_scale_forcing_nudging', 'PA0370', 0, 0, 0, 6, 0 )
1137       ENDIF
1138
1139
1140    END SUBROUTINE nudge_init
1141
1142!------------------------------------------------------------------------------!
1143! Description:
1144! ------------
1145!> @todo Missing subroutine description.
1146!------------------------------------------------------------------------------!
1147    SUBROUTINE calc_tnudge ( time )
1148
1149       IMPLICIT NONE
1150
1151
1152       REAL(wp) ::  dtm         !<
1153       REAL(wp) ::  dtp         !<
1154       REAL(wp) ::  time        !<
1155
1156       INTEGER(iwp) ::  k   !<
1157       INTEGER(iwp) ::  nt  !<
1158
1159       nt = 1
1160       DO WHILE ( time > timenudge(nt) )
1161         nt = nt+1
1162       ENDDO
1163       IF ( time /= timenudge(1) ) THEN
1164         nt = nt-1
1165       ENDIF
1166
1167       dtm = ( time - timenudge(nt) ) / ( timenudge(nt+1) - timenudge(nt) )
1168       dtp = ( timenudge(nt+1) - time ) / ( timenudge(nt+1) - timenudge(nt) )
1169
1170       DO  k = nzb, nzt
1171          tmp_tnudge(k) = MAX( dt_3d, tnudge(k,nt) * dtp + tnudge(k,nt+1) * dtm )
1172       ENDDO
1173
1174    END SUBROUTINE calc_tnudge
1175
1176!------------------------------------------------------------------------------!
1177! Description:
1178! ------------
1179!> Call for all grid points
1180!------------------------------------------------------------------------------!
1181    SUBROUTINE nudge ( time, prog_var )
1182
1183       IMPLICIT NONE
1184
1185       CHARACTER (LEN=*) ::  prog_var  !<
1186
1187       REAL(wp) ::  tmp_tend    !<
1188       REAL(wp) ::  dtm         !<
1189       REAL(wp) ::  dtp         !<
1190       REAL(wp) ::  time        !<
1191
1192       INTEGER(iwp) ::  i  !<
1193       INTEGER(iwp) ::  j  !<
1194       INTEGER(iwp) ::  k  !<
1195       INTEGER(iwp) ::  nt  !<
1196
1197
1198       nt = 1
1199       DO WHILE ( time > timenudge(nt) )
1200         nt = nt+1
1201       ENDDO
1202       IF ( time /= timenudge(1) ) THEN
1203         nt = nt-1
1204       ENDIF
1205
1206       dtm = ( time - timenudge(nt) ) / ( timenudge(nt+1) - timenudge(nt) )
1207       dtp = ( timenudge(nt+1) - time ) / ( timenudge(nt+1) - timenudge(nt) )
1208
1209       SELECT CASE ( prog_var )
1210
1211          CASE ( 'u' )
1212
1213             DO  i = nxl, nxr
1214                DO  j = nys, nyn
1215
1216                   DO  k = nzb+1, nzt
1217
1218                      tmp_tend = - ( hom(k,1,1,0) - ( unudge(k,nt) * dtp +     &
1219                                     unudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
1220
1221                      tend(k,j,i) = tend(k,j,i) + tmp_tend *                   &
1222                                        MERGE( 1.0_wp, 0.0_wp,                 &
1223                                               BTEST( wall_flags_0(k,j,i), 1 ) )
1224
1225                      sums_ls_l(k,6) = sums_ls_l(k,6) + tmp_tend *             &
1226                                     weight_substep(intermediate_timestep_count)
1227                   ENDDO
1228                 
1229                   sums_ls_l(nzt+1,6) = sums_ls_l(nzt,6)
1230 
1231                ENDDO
1232            ENDDO
1233
1234          CASE ( 'v' )
1235
1236             DO  i = nxl, nxr
1237                DO  j = nys, nyn
1238
1239                   DO  k = nzb+1, nzt
1240
1241                      tmp_tend = - ( hom(k,1,2,0) - ( vnudge(k,nt) * dtp +     &
1242                                     vnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
1243
1244                      tend(k,j,i) = tend(k,j,i) + tmp_tend *                   &
1245                                        MERGE( 1.0_wp, 0.0_wp,                 &
1246                                               BTEST( wall_flags_0(k,j,i), 2 ) )
1247
1248                      sums_ls_l(k,7) = sums_ls_l(k,7) + tmp_tend *             &
1249                                     weight_substep(intermediate_timestep_count)
1250                   ENDDO
1251                 
1252                   sums_ls_l(nzt+1,7) = sums_ls_l(nzt,7)
1253
1254                ENDDO
1255            ENDDO
1256
1257          CASE ( 'pt' )
1258
1259             DO  i = nxl, nxr
1260                DO  j = nys, nyn
1261
1262                   DO  k = nzb+1, nzt
1263
1264                      tmp_tend = - ( hom(k,1,4,0) - ( ptnudge(k,nt) * dtp +    &
1265                                     ptnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
1266
1267                      tend(k,j,i) = tend(k,j,i) + tmp_tend *                   &
1268                                        MERGE( 1.0_wp, 0.0_wp,                 &
1269                                               BTEST( wall_flags_0(k,j,i), 0 ) )
1270
1271                      sums_ls_l(k,4) = sums_ls_l(k,4) + tmp_tend *             &
1272                                     weight_substep(intermediate_timestep_count)
1273                   ENDDO
1274
1275                   sums_ls_l(nzt+1,4) = sums_ls_l(nzt,4)
1276
1277                ENDDO
1278            ENDDO
1279
1280          CASE ( 'q' )
1281
1282             DO  i = nxl, nxr
1283                DO  j = nys, nyn
1284
1285                   DO  k = nzb+1, nzt
1286
1287                      tmp_tend = - ( hom(k,1,41,0) - ( qnudge(k,nt) * dtp +    &
1288                                     qnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
1289
1290                      tend(k,j,i) = tend(k,j,i) + tmp_tend *                   &
1291                                        MERGE( 1.0_wp, 0.0_wp,                 &
1292                                               BTEST( wall_flags_0(k,j,i), 0 ) )
1293
1294                      sums_ls_l(k,5) = sums_ls_l(k,5) + tmp_tend *             &
1295                                     weight_substep(intermediate_timestep_count)
1296                   ENDDO
1297                 
1298                   sums_ls_l(nzt+1,5) = sums_ls_l(nzt,5)
1299
1300                ENDDO
1301            ENDDO
1302
1303          CASE DEFAULT
1304             message_string = 'unknown prognostic variable "' // prog_var // '"'
1305             CALL message( 'nudge', 'PA0367', 1, 2, 0, 6, 0 )
1306
1307       END SELECT
1308
1309    END SUBROUTINE nudge
1310
1311
1312!------------------------------------------------------------------------------!
1313! Description:
1314! ------------
1315!> Call for grid point i,j
1316!------------------------------------------------------------------------------!
1317
1318    SUBROUTINE nudge_ij( i, j, time, prog_var )
1319
1320       IMPLICIT NONE
1321
1322
1323       CHARACTER (LEN=*) ::  prog_var  !<
1324
1325       REAL(wp) ::  tmp_tend    !<
1326       REAL(wp) ::  dtm         !<
1327       REAL(wp) ::  dtp         !<
1328       REAL(wp) ::  time        !<
1329
1330       INTEGER(iwp) ::  i  !<
1331       INTEGER(iwp) ::  j  !<
1332       INTEGER(iwp) ::  k  !<
1333       INTEGER(iwp) ::  nt  !<
1334
1335
1336       nt = 1
1337       DO WHILE ( time > timenudge(nt) )
1338         nt = nt+1
1339       ENDDO
1340       IF ( time /= timenudge(1) )  THEN
1341         nt = nt-1
1342       ENDIF
1343
1344       dtm = ( time - timenudge(nt) ) / ( timenudge(nt+1) - timenudge(nt) )
1345       dtp = ( timenudge(nt+1) - time ) / ( timenudge(nt+1) - timenudge(nt) )
1346
1347       SELECT CASE ( prog_var )
1348
1349          CASE ( 'u' )
1350
1351             DO  k = nzb+1, nzt
1352
1353                tmp_tend = - ( hom(k,1,1,0) - ( unudge(k,nt) * dtp +           &
1354                               unudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
1355
1356                tend(k,j,i) = tend(k,j,i) + tmp_tend *                         &
1357                                        MERGE( 1.0_wp, 0.0_wp,                 &
1358                                               BTEST( wall_flags_0(k,j,i), 1 ) )
1359
1360                sums_ls_l(k,6) = sums_ls_l(k,6) + tmp_tend                     &
1361                                 * weight_substep(intermediate_timestep_count)
1362             ENDDO
1363
1364             sums_ls_l(nzt+1,6) = sums_ls_l(nzt,6)
1365
1366          CASE ( 'v' )
1367
1368             DO  k = nzb+1, nzt
1369
1370                tmp_tend = - ( hom(k,1,2,0) - ( vnudge(k,nt) * dtp +           &
1371                               vnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
1372
1373                tend(k,j,i) = tend(k,j,i) + tmp_tend *                         &
1374                                        MERGE( 1.0_wp, 0.0_wp,                 &
1375                                               BTEST( wall_flags_0(k,j,i), 2 ) )
1376
1377                sums_ls_l(k,7) = sums_ls_l(k,7) + tmp_tend                     &
1378                                 * weight_substep(intermediate_timestep_count)
1379             ENDDO
1380
1381             sums_ls_l(nzt+1,7) = sums_ls_l(nzt,7)
1382
1383          CASE ( 'pt' )
1384
1385             DO  k = nzb+1, nzt
1386
1387                tmp_tend = - ( hom(k,1,4,0) - ( ptnudge(k,nt) * dtp +          &
1388                               ptnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
1389
1390                tend(k,j,i) = tend(k,j,i) + tmp_tend *                         &
1391                                        MERGE( 1.0_wp, 0.0_wp,                 &
1392                                               BTEST( wall_flags_0(k,j,i), 0 ) )
1393
1394                sums_ls_l(k,4) = sums_ls_l(k,4) + tmp_tend                     &
1395                                 * weight_substep(intermediate_timestep_count)
1396             ENDDO
1397
1398             sums_ls_l(nzt+1,4) = sums_ls_l(nzt,4)
1399
1400
1401          CASE ( 'q' )
1402
1403             DO  k = nzb+1, nzt
1404
1405                tmp_tend = - ( hom(k,1,41,0) - ( qnudge(k,nt) * dtp +          &
1406                               qnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
1407
1408                tend(k,j,i) = tend(k,j,i) + tmp_tend *                         &
1409                                        MERGE( 1.0_wp, 0.0_wp,                 &
1410                                               BTEST( wall_flags_0(k,j,i), 0 ) )
1411
1412                sums_ls_l(k,5) = sums_ls_l(k,5) + tmp_tend                     &
1413                                 * weight_substep(intermediate_timestep_count)
1414             ENDDO
1415
1416             sums_ls_l(nzt+1,5) = sums_ls_l(nzt,5)
1417
1418          CASE DEFAULT
1419             message_string = 'unknown prognostic variable "' // prog_var // '"'
1420             CALL message( 'nudge', 'PA0367', 1, 2, 0, 6, 0 )
1421
1422       END SELECT
1423
1424
1425    END SUBROUTINE nudge_ij
1426
1427
1428!------------------------------------------------------------------------------!
1429! Description:
1430! ------------
1431!> @todo Missing subroutine description.
1432!------------------------------------------------------------------------------!
1433    SUBROUTINE nudge_ref ( time )
1434
1435       IMPLICIT NONE
1436
1437       INTEGER(iwp) ::  nt                    !<
1438
1439       REAL(wp)             ::  fac           !<
1440       REAL(wp), INTENT(in) ::  time          !<
1441
1442!
1443!--    Interpolation in time of NUDGING_DATA for pt_init and q_init. This is
1444!--    needed for correct upper boundary conditions for pt and q and in case that
1445!      large scale subsidence as well as scalar Rayleigh-damping are used
1446       nt = 1
1447       DO WHILE ( time > time_vert(nt) )
1448          nt = nt + 1
1449       ENDDO
1450       IF ( time /= time_vert(nt) )  THEN
1451        nt = nt - 1
1452       ENDIF
1453
1454       fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
1455
1456       pt_init = ptnudge(:,nt) + fac * ( ptnudge(:,nt+1) - ptnudge(:,nt) )
1457       q_init  = qnudge(:,nt) + fac * ( qnudge(:,nt+1) - qnudge(:,nt) )
1458       u_init  = unudge(:,nt) + fac * ( unudge(:,nt+1) - unudge(:,nt) )
1459       v_init  = vnudge(:,nt) + fac * ( vnudge(:,nt+1) - vnudge(:,nt) )
1460
1461    END SUBROUTINE nudge_ref
1462
1463
1464 END MODULE lsf_nudging_mod
Note: See TracBrowser for help on using the repository browser.