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

Last change on this file since 2320 was 2320, checked in by suehring, 7 years ago

large-scale forcing and nudging modularized

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