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

Last change on this file since 4867 was 4828, checked in by Giersch, 4 years ago

Copyright updated to year 2021, interface pmc_sort removed to accelarate the nesting code

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