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

Last change on this file since 2350 was 2342, checked in by boeske, 7 years ago

fixed check if surface forcing data is available until end of simulation

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