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

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