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

Last change on this file since 4428 was 4360, checked in by suehring, 4 years ago

Bugfix in output of time-averaged plant-canopy quanities; Output of plant-canopy data only where tall canopy is defined; land-surface model: fix wrong location strings; tests: update urban test case; all source code files: copyright update

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