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

Last change on this file since 3582 was 3428, checked in by gronemeier, 5 years ago

Rename td_XXX_lpt to td_XXX_thetal

  • Property svn:keywords set to Id
File size: 55.4 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 3428 2018-10-25 12:32:05Z suehring $
27! Rename td_XXX_lpt to td_XXX_thetal
28!
29! 3347 2018-10-15 14:21:08Z suehring
30! Offline nesting parts are moved to an independent module nesting_offl_mod
31!
32! 3341 2018-10-15 10:31:27Z suehring
33! ocean renamed ocean_mode
34!
35! 3241 2018-09-12 15:02:00Z raasch
36! unused variables removed
37!
38! 3183 2018-07-27 14:25:55Z suehring
39! * Adjustment to new Inifor version:
40!   - No vertical interpolation/extrapolation of lateral boundary data required
41!     any more (Inifor can treat grid stretching now
42!   - Revise initialization in case of COSMO forcing
43! * Rename variables and subroutines for offline nesting
44!
45! 3182 2018-07-27 13:36:03Z suehring
46! Error messages revised
47!
48! 3045 2018-05-28 07:55:41Z Giersch
49! Error messages revised
50!
51! 3026 2018-05-22 10:30:53Z schwenkel
52! Changed the name specific humidity to mixing ratio, since we are computing
53! mixing ratios.
54!
55! 2970 2018-04-13 15:09:23Z suehring
56! Bugfix in old large-scale forcing mode
57!
58! 2938 2018-03-27 15:52:42Z suehring
59! Further improvements for nesting in larger-scale model
60!
61! 2863 2018-03-08 11:36:25Z suehring
62! Corrected "Former revisions" section
63!
64! 2696 2017-12-14 17:12:51Z kanani
65! Change in file header (GPL part)
66! Forcing with larger-scale models implemented (MS)
67!
68! 2342 2017-08-08 11:00:43Z boeske
69! fixed check if surface forcing data is available until end of simulation
70!
71! 2320 2017-07-21 12:47:43Z suehring
72! initial revision
73!
74! Description:
75! ------------
76!> Calculates large scale forcings (geostrophic wind and subsidence velocity) as
77!> well as surfaces fluxes dependent on time given in an external file (LSF_DATA).
78!> Moreover, module contains nudging routines, where u, v, pt and q are nudged
79!> to given profiles on a relaxation timescale tnudge.
80!> Profiles are read in from NUDGING_DATA.
81!> Code is based on Neggers et al. (2012) and also in parts on DALES and UCLA-LES.
82!> @todo: Revise reading of ASCII-files
83!> @todo: Remove unused variables and control flags
84!> @todo: Revise large-scale facing of surface variables
85!> @todo: Revise control flags lsf_exception, lsf_surf, lsf_vert, etc.
86!--------------------------------------------------------------------------------!
87 MODULE lsf_nudging_mod
88
89    USE arrays_3d,                                                             &
90        ONLY:  dzw, e, diss, heatflux_input_conversion, pt, pt_init, q,        &
91               q_init, s, tend, u, u_init, ug, v, v_init, vg, w, w_subs,       &
92               waterflux_input_conversion, zu, zw                 
93
94    USE control_parameters,                                                    &
95        ONLY:  bc_lr, bc_ns, bc_pt_b, bc_q_b, constant_diffusion,              &
96               constant_heatflux, constant_waterflux,                          &
97               data_output_pr, dt_3d, end_time,                                &
98               humidity, initializing_actions, intermediate_timestep_count,    &
99               ibc_pt_b, ibc_q_b,                                              &
100               large_scale_forcing, large_scale_subsidence, lsf_surf, lsf_vert,&
101               lsf_exception, message_string, neutral,                         &
102               nudging, passive_scalar, pt_surface, ocean_mode, q_surface,     &
103               surface_heatflux, surface_pressure, surface_waterflux,          &
104               topography, use_subsidence_tendencies
105               
106    USE cpulog,                                                                &
107        ONLY:  cpu_log, log_point
108
109    USE grid_variables
110
111    USE indices,                                                               &
112        ONLY:  nbgp, ngp_sums_ls, nx, nxl, nxlg, nxlu, nxr, nxrg, ny, nys,     &
113               nysv, nysg, nyn, nyng, nzb, nz, nzt, wall_flags_0
114
115    USE kinds
116
117    USE pegrid
118
119    USE surface_mod,                                                           &
120        ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
121
122    USE statistics,                                                            &
123        ONLY:  hom, statistic_regions, sums_ls_l, weight_substep
124
125    INTEGER(iwp) ::  nlsf = 1000                       !< maximum number of profiles in LSF_DATA (large scale forcing)
126    INTEGER(iwp) ::  ntnudge = 1000                    !< maximum number of profiles in NUDGING_DATA (nudging)
127
128    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ptnudge     !< vertical profile of pot. temperature interpolated to vertical grid (nudging)
129    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  qnudge      !< vertical profile of water vapor mixing ratio interpolated to vertical grid (nudging)
130    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tnudge      !< vertical profile of nudging time scale interpolated to vertical grid (nudging) 
131    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  td_lsa_lpt  !< temperature tendency due to large scale advection (large scale forcing)
132    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  td_lsa_q    !< water vapor mixing ratio tendency due to large scale advection (large scale forcing)
133    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  td_sub_lpt  !< temperature tendency due to subsidence/ascent (large scale forcing)
134    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  td_sub_q    !< water vapor mixing ratio tendency due to subsidence/ascent (large scale forcing)
135    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ug_vert     !< vertical profile of geostrophic wind component in x-direction interpolated to vertical grid (large scale forcing)
136    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  unudge      !< vertical profile of wind component in x-direction interpolated to vertical grid (nudging)
137    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  vnudge      !< vertical profile of wind component in y-direction interpolated to vertical grid (nudging)
138    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  vg_vert     !< vertical profile of geostrophic wind component in y-direction interpolated to vertical grid (large scale forcing)
139    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  wnudge      !< vertical profile of subsidence/ascent velocity interpolated to vertical grid (nudging) ???
140    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  wsubs_vert  !< vertical profile of wind component in z-direction interpolated to vertical grid (nudging) ???
141
142    REAL(wp), DIMENSION(:), ALLOCATABLE ::  shf_surf      !< time-dependent surface sensible heat flux (large scale forcing)
143    REAL(wp), DIMENSION(:), ALLOCATABLE ::  timenudge     !< times at which vertical profiles are defined in NUDGING_DATA (nudging)
144    REAL(wp), DIMENSION(:), ALLOCATABLE ::  time_surf     !< times at which surface values/fluxes are defined in LSF_DATA (large scale forcing)
145    REAL(wp), DIMENSION(:), ALLOCATABLE ::  time_vert     !< times at which vertical profiles are defined in LSF_DATA (large scale forcing)
146
147    REAL(wp), DIMENSION(:), ALLOCATABLE ::  tmp_tnudge    !< current nudging time scale
148
149    REAL(wp), DIMENSION(:), ALLOCATABLE ::  p_surf        !< time-dependent surface pressure (large scale forcing)
150    REAL(wp), DIMENSION(:), ALLOCATABLE ::  pt_surf       !< time-dependent surface temperature (large scale forcing)
151    REAL(wp), DIMENSION(:), ALLOCATABLE ::  qsws_surf     !< time-dependent surface latent heat flux (large scale forcing)
152    REAL(wp), DIMENSION(:), ALLOCATABLE ::  q_surf        !< time-dependent surface water vapor mixing ratio (large scale forcing)
153
154    SAVE
155    PRIVATE
156!
157!-- Public subroutines
158    PUBLIC calc_tnudge, ls_forcing_surf, ls_forcing_vert, ls_advec, lsf_init,  &
159           lsf_nudging_check_parameters, nudge_init,                           &
160           lsf_nudging_check_data_output_pr, lsf_nudging_header,               &
161           nudge, nudge_ref
162           
163!
164!-- Public variables
165    PUBLIC qsws_surf, shf_surf, td_lsa_lpt, td_lsa_q, td_sub_lpt,              &
166           td_sub_q, time_vert
167
168
169    INTERFACE ls_advec
170       MODULE PROCEDURE ls_advec
171       MODULE PROCEDURE ls_advec_ij
172    END INTERFACE ls_advec
173
174    INTERFACE nudge
175       MODULE PROCEDURE nudge
176       MODULE PROCEDURE nudge_ij
177    END INTERFACE nudge
178
179 CONTAINS
180
181
182!------------------------------------------------------------------------------!
183! Description:
184! ------------
185!> @todo Missing subroutine description.
186!------------------------------------------------------------------------------!
187    SUBROUTINE lsf_nudging_check_parameters
188
189       IMPLICIT NONE
190!
191!--    Check nudging and large scale forcing from external file
192       IF ( nudging  .AND.  (  .NOT.  large_scale_forcing ) )  THEN
193          message_string = 'Nudging requires large_scale_forcing = .T.. &'//   &
194                        'Surface fluxes and geostrophic wind should be &'//    &
195                        'prescribed in file LSF_DATA'
196          CALL message( 'check_parameters', 'PA0374', 1, 2, 0, 6, 0 )
197       ENDIF
198
199       IF ( large_scale_forcing  .AND.  ( bc_lr /= 'cyclic'  .OR.              &
200                                          bc_ns /= 'cyclic' ) )  THEN
201          message_string = 'Non-cyclic lateral boundaries do not allow for &'//&
202                        'the usage of large scale forcing from external file.'
203          CALL message( 'check_parameters', 'PA0375', 1, 2, 0, 6, 0 )
204       ENDIF
205
206       IF ( large_scale_forcing  .AND.  (  .NOT.  humidity ) )  THEN
207          message_string = 'The usage of large scale forcing from external &'//&
208                        'file LSF_DATA requires humidity = .T..'
209          CALL message( 'check_parameters', 'PA0376', 1, 2, 0, 6, 0 )
210       ENDIF
211
212       IF ( large_scale_forcing  .AND.  passive_scalar )  THEN
213          message_string = 'The usage of large scale forcing from external &'// &
214                        'file LSF_DATA is not implemented for passive scalars'
215          CALL message( 'check_parameters', 'PA0440', 1, 2, 0, 6, 0 )
216       ENDIF
217
218       IF ( large_scale_forcing  .AND.  topography /= 'flat'                   &
219                              .AND.  .NOT.  lsf_exception )  THEN
220          message_string = 'The usage of large scale forcing from external &'//&
221                        'file LSF_DATA is not implemented for non-flat topography'
222          CALL message( 'check_parameters', 'PA0377', 1, 2, 0, 6, 0 )
223       ENDIF
224
225       IF ( large_scale_forcing  .AND.  ocean_mode )  THEN
226          message_string = 'The usage of large scale forcing from external &'//&
227                        'file LSF_DATA is not implemented for ocean mode'
228          CALL message( 'check_parameters', 'PA0378', 1, 2, 0, 6, 0 )
229       ENDIF
230
231    END SUBROUTINE lsf_nudging_check_parameters
232
233!------------------------------------------------------------------------------!
234! Description:
235! ------------
236!> Check data output of profiles for land surface model
237!------------------------------------------------------------------------------!
238    SUBROUTINE lsf_nudging_check_data_output_pr( variable, var_count, unit,    &
239                                                 dopr_unit )
240 
241       USE profil_parameter
242
243       IMPLICIT NONE
244   
245       CHARACTER (LEN=*) ::  unit      !<
246       CHARACTER (LEN=*) ::  variable  !<
247       CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
248 
249       INTEGER(iwp) ::  var_count     !<
250
251       SELECT CASE ( TRIM( variable ) )
252       
253
254          CASE ( 'td_lsa_thetal' )
255             IF (  .NOT.  large_scale_forcing )  THEN
256                message_string = 'data_output_pr = ' //                        &
257                                 TRIM( data_output_pr(var_count) ) //          &
258                                 ' is not implemented for ' //                 &
259                                 'large_scale_forcing = .FALSE.'
260                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393',    &
261                               1, 2, 0, 6, 0 )
262             ELSE
263                dopr_index(var_count) = 81
264                dopr_unit             = 'K/s'
265                unit                  = 'K/s'
266                hom(:,2,81,:) = SPREAD( zu, 2, statistic_regions+1 )
267             ENDIF
268
269          CASE ( 'td_lsa_q' )
270             IF (  .NOT.  large_scale_forcing )  THEN
271                message_string = 'data_output_pr = ' //                        &
272                                 TRIM( data_output_pr(var_count) ) //          &
273                                 ' is not implemented for ' //                 &
274                                 'large_scale_forcing = .FALSE.'
275                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393',    &
276                               1, 2, 0, 6, 0 )
277             ELSE
278                dopr_index(var_count) = 82
279                dopr_unit             = 'kg/kgs'
280                unit                  = 'kg/kgs'
281                hom(:,2,82,:) = SPREAD( zu, 2, statistic_regions+1 )
282             ENDIF
283          CASE ( 'td_sub_thetal' )
284             IF (  .NOT.  large_scale_forcing )  THEN
285                message_string = 'data_output_pr = ' //                        &
286                                 TRIM( data_output_pr(var_count) ) //          &
287                                 ' is not implemented for ' //                 &
288                                 'large_scale_forcing = .FALSE.'
289                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393',    &
290                               1, 2, 0, 6, 0 )
291             ELSE
292                dopr_index(var_count) = 83
293                dopr_unit             = 'K/s'
294                unit                  = 'K/s'
295                hom(:,2,83,:) = SPREAD( zu, 2, statistic_regions+1 )
296             ENDIF
297
298          CASE ( 'td_sub_q' )
299             IF (  .NOT.  large_scale_forcing )  THEN
300                message_string = 'data_output_pr = ' //                        &
301                                 TRIM( data_output_pr(var_count) ) //          &
302                                 ' is not implemented for ' //                 &
303                                 'large_scale_forcing = .FALSE.'
304                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393',    &
305                               1, 2, 0, 6, 0 )
306             ELSE
307                dopr_index(var_count) = 84
308                dopr_unit             = 'kg/kgs'
309                unit                  = 'kg/kgs'
310                hom(:,2,84,:) = SPREAD( zu, 2, statistic_regions+1 )
311             ENDIF
312
313          CASE ( 'td_nud_thetal' )
314             IF (  .NOT.  nudging )  THEN
315                message_string = 'data_output_pr = ' //                        &
316                                 TRIM( data_output_pr(var_count) ) //          &
317                                 ' is not implemented for ' //                 &
318                                 'nudging = .FALSE.'
319                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0394',    &
320                               1, 2, 0, 6, 0 )
321             ELSE
322                dopr_index(var_count) = 85
323                dopr_unit             = 'K/s'
324                unit                  = 'K/s'
325                hom(:,2,85,:) = SPREAD( zu, 2, statistic_regions+1 )
326             ENDIF
327
328          CASE ( 'td_nud_q' )
329             IF (  .NOT.  nudging )  THEN
330                message_string = 'data_output_pr = ' //                        &
331                                 TRIM( data_output_pr(var_count) ) //          &
332                                 ' is not implemented for ' //                 &
333                                 'nudging = .FALSE.'
334                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0394',    &
335                               1, 2, 0, 6, 0 )
336             ELSE
337                dopr_index(var_count) = 86
338                dopr_unit             = 'kg/kgs'
339                unit                  = 'kg/kgs'
340                hom(:,2,86,:) = SPREAD( zu, 2, statistic_regions+1 )
341             ENDIF
342
343          CASE ( 'td_nud_u' )
344             IF (  .NOT.  nudging )  THEN
345                message_string = 'data_output_pr = ' //                        &
346                                 TRIM( data_output_pr(var_count) ) //          &
347                                 ' is not implemented for ' //                 &
348                                 'nudging = .FALSE.'
349                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0394',    &
350                               1, 2, 0, 6, 0 )
351             ELSE
352                dopr_index(var_count) = 87
353                dopr_unit             = 'm/s2'
354                unit                  = 'm/s2'
355                hom(:,2,87,:) = SPREAD( zu, 2, statistic_regions+1 )
356             ENDIF
357
358          CASE ( 'td_nud_v' )
359             IF (  .NOT.  nudging )  THEN
360                message_string = 'data_output_pr = ' //                        &
361                                 TRIM( data_output_pr(var_count) ) //          &
362                                 ' is not implemented for ' //                 &
363                                 'nudging = .FALSE.'
364                CALL message( 'lsf_nudging_check_data_output_pr', 'PA0394',    &
365                               1, 2, 0, 6, 0 )
366             ELSE
367                dopr_index(var_count) = 88
368                dopr_unit             = 'm/s2'
369                unit                  = 'm/s2'
370                hom(:,2,88,:) = SPREAD( zu, 2, statistic_regions+1 )
371             ENDIF
372
373
374          CASE DEFAULT
375             unit = 'illegal'
376   
377       END SELECT
378
379    END SUBROUTINE lsf_nudging_check_data_output_pr
380
381!------------------------------------------------------------------------------!
382! Description:
383! ------------
384!> @todo Missing subroutine description.
385!------------------------------------------------------------------------------!
386    SUBROUTINE lsf_nudging_header ( io )
387
388       IMPLICIT NONE
389
390       INTEGER(iwp), INTENT(IN) ::  io !< Unit of the output file
391
392       WRITE ( io, 1 )
393       IF ( large_scale_forcing )  THEN
394          WRITE ( io, 3 )
395          WRITE ( io, 4 )
396
397          IF ( large_scale_subsidence )  THEN
398             IF ( .NOT. use_subsidence_tendencies )  THEN
399                WRITE ( io, 5 )
400             ELSE
401                WRITE ( io, 6 )
402             ENDIF
403          ENDIF
404
405          IF ( bc_pt_b == 'dirichlet' )  THEN
406             WRITE ( io, 12 )
407          ELSEIF ( bc_pt_b == 'neumann' )  THEN
408             WRITE ( io, 13 )
409          ENDIF
410
411          IF ( bc_q_b == 'dirichlet' )  THEN
412             WRITE ( io, 14 )
413          ELSEIF ( bc_q_b == 'neumann' )  THEN
414             WRITE ( io, 15 )
415          ENDIF
416
417          WRITE ( io, 7 )
418          IF ( nudging )  THEN
419             WRITE ( io, 10 )
420          ENDIF
421       ELSE
422          WRITE ( io, 2 )
423          WRITE ( io, 11 )
424       ENDIF
425       IF ( large_scale_subsidence )  THEN
426          WRITE ( io, 8 )
427          WRITE ( io, 9 )
428       ENDIF
429
430
4311 FORMAT (//' Large scale forcing and nudging:'/ &
432              ' -------------------------------'/)
4332 FORMAT (' --> No large scale forcing from external is used (default) ')
4343 FORMAT (' --> Large scale forcing from external file LSF_DATA is used: ')
4354 FORMAT ('     - large scale advection tendencies ')
4365 FORMAT ('     - large scale subsidence velocity w_subs ')
4376 FORMAT ('     - large scale subsidence tendencies ')
4387 FORMAT ('     - and geostrophic wind components ug and vg')
4398 FORMAT (' --> Large-scale vertical motion is used in the ', &
440                  'prognostic equation(s) for')
4419 FORMAT ('     the scalar(s) only')
44210 FORMAT (' --> Nudging is used')
44311 FORMAT (' --> No nudging is used (default) ')
44412 FORMAT ('     - prescribed surface values for temperature')
44513 FORMAT ('     - prescribed surface fluxes for temperature')
44614 FORMAT ('     - prescribed surface values for humidity')
44715 FORMAT ('     - prescribed surface fluxes for humidity')
448
449    END SUBROUTINE lsf_nudging_header 
450
451!------------------------------------------------------------------------------!
452! Description:
453! ------------
454!> @todo Missing subroutine description.
455!------------------------------------------------------------------------------!
456    SUBROUTINE lsf_init
457
458       USE control_parameters,                                                 &
459           ONLY:  bc_lr_cyc, bc_ns_cyc
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.