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

Last change on this file since 4180 was 4180, checked in by scharf, 2 years ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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