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

Last change on this file since 3395 was 3347, checked in by suehring, 5 years ago

Offline nesting revised and separated from large_scale_forcing_mod; Time-dependent synthetic turbulence generator; bugfixes in USM/LSM radiation model and init_pegrid

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