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

Last change on this file since 4332 was 4329, checked in by motisi, 4 years ago

Renamed wall_flags_0 to wall_flags_static_0

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