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

Last change on this file since 2350 was 2339, checked in by gronemeier, 7 years ago

corrected timestamp in header

  • Property svn:keywords set to Id
File size: 37.1 KB
Line 
1!> @file model_1d_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: model_1d_mod.f90 2339 2017-08-07 13:55:26Z kanani $
27! corrected timestamp in header
28!
29! 2338 2017-08-07 12:15:38Z gronemeier
30! renamed init_1d_model to model_1d_mod and and formatted it as a module;
31! reformatted output of profiles
32!
33! 2337 2017-08-07 08:59:53Z gronemeier
34! revised calculation of mixing length
35! removed rounding of time step
36! corrected calculation of virtual potential temperature
37!
38! 2334 2017-08-04 11:57:04Z gronemeier
39! set c_m = 0.4 according to Detering and Etling (1985)
40!
41! 2299 2017-06-29 10:14:38Z maronga
42! Removed german text
43!
44! 2101 2017-01-05 16:42:31Z suehring
45!
46! 2059 2016-11-10 14:20:40Z maronga
47! Corrected min/max values of Rif.
48!
49! 2000 2016-08-20 18:09:15Z knoop
50! Forced header and separation lines into 80 columns
51!
52! 1960 2016-07-12 16:34:24Z suehring
53! Remove passive_scalar from IF-statements, as 1D-scalar profile is effectively
54! not used.
55! Formatting adjustment
56!
57! 1808 2016-04-05 19:44:00Z raasch
58! routine local_flush replaced by FORTRAN statement
59!
60! 1709 2015-11-04 14:47:01Z maronga
61! Set initial time step to 10 s to avoid instability of the 1d model for small
62! grid spacings
63!
64! 1697 2015-10-28 17:14:10Z raasch
65! small E- and F-FORMAT changes to avoid informative compiler messages about
66! insufficient field width
67!
68! 1691 2015-10-26 16:17:44Z maronga
69! Renamed prandtl_layer to constant_flux_layer. rif is replaced by ol and zeta.
70!
71! 1682 2015-10-07 23:56:08Z knoop
72! Code annotations made doxygen readable
73!
74! 1353 2014-04-08 15:21:23Z heinze
75! REAL constants provided with KIND-attribute
76!
77! 1346 2014-03-27 13:18:20Z heinze
78! Bugfix: REAL constants provided with KIND-attribute especially in call of
79! intrinsic function like MAX, MIN, SIGN
80!
81! 1322 2014-03-20 16:38:49Z raasch
82! REAL functions provided with KIND-attribute
83!
84! 1320 2014-03-20 08:40:49Z raasch
85! ONLY-attribute added to USE-statements,
86! kind-parameters added to all INTEGER and REAL declaration statements,
87! kinds are defined in new module kinds,
88! revision history before 2012 removed,
89! comment fields (!:) to be used for variable explanations added to
90! all variable declaration statements
91!
92! 1036 2012-10-22 13:43:42Z raasch
93! code put under GPL (PALM 3.9)
94!
95! 1015 2012-09-27 09:23:24Z raasch
96! adjustment of mixing length to the Prandtl mixing length at first grid point
97! above ground removed
98!
99! 1001 2012-09-13 14:08:46Z raasch
100! all actions concerning leapfrog scheme removed
101!
102! 996 2012-09-07 10:41:47Z raasch
103! little reformatting
104!
105! 978 2012-08-09 08:28:32Z fricke
106! roughness length for scalar quantities z0h1d added
107!
108! Revision 1.1  1998/03/09 16:22:10  raasch
109! Initial revision
110!
111!
112! Description:
113! ------------
114!> 1D-model to initialize the 3D-arrays.
115!> The temperature profile is set as steady and a corresponding steady solution
116!> of the wind profile is being computed.
117!> All subroutines required can be found within this file.
118!>
119!> @todo harmonize code with new surface_layer_fluxes module
120!> @bug 1D model crashes when using small grid spacings in the order of 1 m
121!------------------------------------------------------------------------------!
122 MODULE model_1d_mod
123
124    USE arrays_3d,                                                             &
125        ONLY:  dd2zu, ddzu, ddzw, dzu, l_grid, pt_init, q_init, ug, u_init,    &
126               vg, v_init, zu
127   
128    USE control_parameters,                                                    &
129        ONLY:  constant_diffusion, constant_flux_layer, dissipation_1d, f, g,  &
130               humidity, ibc_e_b, intermediate_timestep_count,                 &
131               intermediate_timestep_count_max, kappa, km_constant,            &
132               message_string, mixing_length_1d, prandtl_number,               &
133               roughness_length, run_description_header, simulated_time_chr, timestep_scheme, tsc, z0h_factor
134
135    USE indices,                                                               &
136        ONLY:  nzb, nzb_diff, nzt
137   
138    USE kinds
139
140    USE pegrid
141       
142
143    IMPLICIT NONE
144
145
146    INTEGER(iwp) ::  current_timestep_number_1d = 0  !< current timestep number (1d-model)
147    INTEGER(iwp) ::  damp_level_ind_1d               !< lower grid index of damping layer (1d-model) !#
148
149    LOGICAL ::  run_control_header_1d = .FALSE.  !< flag for output of run control header (1d-model)
150    LOGICAL ::  stop_dt_1d = .FALSE.             !< termination flag, used in case of too small timestep (1d-model)
151
152    REAL(wp) ::  c_m = 0.4_wp                  !< model constant, 0.4 according to Detering and Etling (1985)
153    REAL(wp) ::  damp_level_1d = -1.0_wp       !< namelist parameter    !#
154    REAL(wp) ::  dt_1d = 60.0_wp               !< dynamic timestep (1d-model)
155    REAL(wp) ::  dt_max_1d = 300.0_wp          !< timestep limit (1d-model)
156    REAL(wp) ::  dt_pr_1d = 9999999.9_wp       !< namelist parameter     !#
157    REAL(wp) ::  dt_run_control_1d = 60.0_wp   !< namelist parameter     !#
158    REAL(wp) ::  end_time_1d = 864000.0_wp     !< namelist parameter     !#
159    REAL(wp) ::  qs1d                          !< characteristic humidity scale (1d-model)
160    REAL(wp) ::  simulated_time_1d = 0.0_wp    !< updated simulated time (1d-model)
161    REAL(wp) ::  time_pr_1d = 0.0_wp           !< updated simulated time for profile output (1d-model)
162    REAL(wp) ::  time_run_control_1d = 0.0_wp  !< updated simulated time for run-control output (1d-model)
163    REAL(wp) ::  ts1d                          !< characteristic temperature scale (1d-model)
164    REAL(wp) ::  us1d                          !< friction velocity (1d-model)                 !#
165    REAL(wp) ::  usws1d                        !< u-component of the momentum flux (1d-model)  !#
166    REAL(wp) ::  vsws1d                        !< v-component of the momentum flux (1d-model)  !#
167    REAL(wp) ::  z01d                          !< roughness length for momentum (1d-model)
168    REAL(wp) ::  z0h1d                         !< roughness length for scalars (1d-model)
169
170
171    REAL(wp), DIMENSION(:), ALLOCATABLE ::  e1d      !< tke (1d-model)                         !#
172    REAL(wp), DIMENSION(:), ALLOCATABLE ::  e1d_p    !< prognostic value of tke (1d-model)
173    REAL(wp), DIMENSION(:), ALLOCATABLE ::  kh1d     !< turbulent diffusion coefficient for heat (1d-model)    !#
174    REAL(wp), DIMENSION(:), ALLOCATABLE ::  km1d     !< turbulent diffusion coefficient for momentum (1d-model)!#
175    REAL(wp), DIMENSION(:), ALLOCATABLE ::  l_black  !< mixing length Blackadar (1d-model)
176    REAL(wp), DIMENSION(:), ALLOCATABLE ::  l1d      !< mixing length for turbulent diffusion coefficients (1d-model) !#
177    REAL(wp), DIMENSION(:), ALLOCATABLE ::  l1d_diss !< mixing length for dissipation (1d-model)
178    REAL(wp), DIMENSION(:), ALLOCATABLE ::  rif1d    !< Richardson flux number (1d-model)      !#
179    REAL(wp), DIMENSION(:), ALLOCATABLE ::  te_e     !< tendency of e (1d-model)
180    REAL(wp), DIMENSION(:), ALLOCATABLE ::  te_em    !< weighted tendency of e for previous sub-timestep (1d-model)
181    REAL(wp), DIMENSION(:), ALLOCATABLE ::  te_u     !< tendency of u (1d-model)
182    REAL(wp), DIMENSION(:), ALLOCATABLE ::  te_um    !< weighted tendency of u for previous sub-timestep (1d-model)
183    REAL(wp), DIMENSION(:), ALLOCATABLE ::  te_v     !< tendency of v (1d-model)
184    REAL(wp), DIMENSION(:), ALLOCATABLE ::  te_vm    !< weighted tendency of v for previous sub-timestep (1d-model)
185    REAL(wp), DIMENSION(:), ALLOCATABLE ::  u1d      !< u-velocity component (1d-model)        !#
186    REAL(wp), DIMENSION(:), ALLOCATABLE ::  u1d_p    !< prognostic value of u-velocity component (1d-model)
187    REAL(wp), DIMENSION(:), ALLOCATABLE ::  v1d      !< v-velocity component (1d-model)        !#
188    REAL(wp), DIMENSION(:), ALLOCATABLE ::  v1d_p    !< prognostic value of v-velocity component (1d-model)
189
190!
191!-- Initialize 1D model
192    INTERFACE init_1d_model
193       MODULE PROCEDURE init_1d_model
194    END INTERFACE init_1d_model
195
196!
197!-- Print profiles
198    INTERFACE print_1d_model
199       MODULE PROCEDURE print_1d_model
200    END INTERFACE print_1d_model
201
202!
203!-- Print run control information
204    INTERFACE run_control_1d
205       MODULE PROCEDURE run_control_1d
206    END INTERFACE run_control_1d
207
208!
209!-- Main procedure
210    INTERFACE time_integration_1d
211       MODULE PROCEDURE time_integration_1d
212    END INTERFACE time_integration_1d
213
214!
215!-- Calculate time step
216    INTERFACE timestep_1d
217       MODULE PROCEDURE timestep_1d
218    END INTERFACE timestep_1d
219
220    SAVE
221
222    PRIVATE
223!
224!-- Public interfaces
225    PUBLIC  init_1d_model
226
227!
228!-- Public variables
229    PUBLIC  damp_level_1d, damp_level_ind_1d, dt_pr_1d, dt_run_control_1d,     &
230            e1d, end_time_1d, kh1d, km1d, l1d, rif1d, u1d, us1d, usws1d, v1d,  &
231            vsws1d
232
233
234    CONTAINS
235
236 SUBROUTINE init_1d_model
237 
238    IMPLICIT NONE
239
240    CHARACTER (LEN=9) ::  time_to_string  !<
241   
242    INTEGER(iwp) ::  k  !<
243   
244    REAL(wp) ::  lambda !<
245
246!
247!-- Allocate required 1D-arrays
248    ALLOCATE( e1d(nzb:nzt+1), e1d_p(nzb:nzt+1), kh1d(nzb:nzt+1),               &
249              km1d(nzb:nzt+1), l_black(nzb:nzt+1), l1d(nzb:nzt+1),             &
250              l1d_diss(nzb:nzt+1), rif1d(nzb:nzt+1), te_e(nzb:nzt+1),          &
251              te_em(nzb:nzt+1), te_u(nzb:nzt+1), te_um(nzb:nzt+1),             &
252              te_v(nzb:nzt+1), te_vm(nzb:nzt+1), u1d(nzb:nzt+1),               &
253              u1d_p(nzb:nzt+1),  v1d(nzb:nzt+1), v1d_p(nzb:nzt+1) )
254
255!
256!-- Initialize arrays
257    IF ( constant_diffusion )  THEN
258       km1d = km_constant
259       kh1d = km_constant / prandtl_number
260    ELSE
261       e1d = 0.0_wp; e1d_p = 0.0_wp
262       kh1d = 0.0_wp; km1d = 0.0_wp
263       rif1d = 0.0_wp
264!
265!--    Compute the mixing length
266       l_black(nzb) = 0.0_wp
267
268       IF ( TRIM( mixing_length_1d ) == 'blackadar' )  THEN
269!
270!--       Blackadar mixing length
271          IF ( f /= 0.0_wp )  THEN
272             lambda = 2.7E-4_wp * SQRT( ug(nzt+1)**2 + vg(nzt+1)**2 ) /        &
273                               ABS( f ) + 1E-10_wp
274          ELSE
275             lambda = 30.0_wp
276          ENDIF
277
278          DO  k = nzb+1, nzt+1
279             l_black(k) = kappa * zu(k) / ( 1.0_wp + kappa * zu(k) / lambda )
280          ENDDO
281
282       ELSEIF ( TRIM( mixing_length_1d ) == 'as_in_3d_model' )  THEN
283!
284!--       Use the same mixing length as in 3D model
285          l_black(1:nzt) = l_grid
286          l_black(nzt+1) = l_black(nzt)
287
288       ENDIF
289    ENDIF
290    l1d      = l_black
291    l1d_diss = l_black
292    u1d      = u_init
293    u1d_p    = u_init
294    v1d      = v_init
295    v1d_p    = v_init
296
297!
298!-- Set initial horizontal velocities at the lowest grid levels to a very small
299!-- value in order to avoid too small time steps caused by the diffusion limit
300!-- in the initial phase of a run (at k=1, dz/2 occurs in the limiting formula!)
301    u1d(0:1)   = 0.1_wp
302    u1d_p(0:1) = 0.1_wp
303    v1d(0:1)   = 0.1_wp
304    v1d_p(0:1) = 0.1_wp
305
306!
307!-- For u*, theta* and the momentum fluxes plausible values are set
308    IF ( constant_flux_layer )  THEN
309       us1d = 0.1_wp   ! without initial friction the flow would not change
310    ELSE
311       e1d(nzb+1)  = 1.0_wp
312       km1d(nzb+1) = 1.0_wp
313       us1d = 0.0_wp
314    ENDIF
315    ts1d = 0.0_wp
316    usws1d = 0.0_wp
317    vsws1d = 0.0_wp
318    z01d  = roughness_length
319    z0h1d = z0h_factor * z01d 
320    IF ( humidity )  qs1d = 0.0_wp
321
322!
323!-- Tendencies must be preset in order to avoid runtime errors within the
324!-- first Runge-Kutta step
325    te_em = 0.0_wp
326    te_um = 0.0_wp
327    te_vm = 0.0_wp
328
329!
330!-- Set model constant
331    IF ( dissipation_1d == 'as_in_3d_model' )  c_m = 0.1_wp
332
333!
334!-- Set start time in hh:mm:ss - format
335    simulated_time_chr = time_to_string( simulated_time_1d )
336
337!
338!-- Integrate the 1D-model equations using the Runge-Kutta scheme
339    CALL time_integration_1d
340
341
342 END SUBROUTINE init_1d_model
343
344
345
346!------------------------------------------------------------------------------!
347! Description:
348! ------------
349!> Runge-Kutta time differencing scheme for the 1D-model.
350!------------------------------------------------------------------------------!
351 
352 SUBROUTINE time_integration_1d
353
354    IMPLICIT NONE
355
356    CHARACTER (LEN=9) ::  time_to_string  !<
357   
358    INTEGER(iwp) ::  k  !< loop index
359   
360    REAL(wp) ::  a            !< auxiliary variable
361    REAL(wp) ::  b            !< auxiliary variable
362    REAL(wp) ::  dissipation  !< dissipation of TKE
363    REAL(wp) ::  dpt_dz       !< vertical temperature gradient
364    REAL(wp) ::  flux         !< vertical temperature gradient
365    REAL(wp) ::  kmzm         !< Km(z-dz/2)
366    REAL(wp) ::  kmzp         !< Km(z+dz/2)
367    REAL(wp) ::  l_stable     !< mixing length for stable case
368    REAL(wp) ::  pt_0         !< reference temperature
369    REAL(wp) ::  uv_total     !< horizontal wind speed
370
371!
372!-- Determine the time step at the start of a 1D-simulation and
373!-- determine and printout quantities used for run control
374    dt_1d = 10.0_wp
375    CALL run_control_1d
376
377!
378!-- Start of time loop
379    DO  WHILE ( simulated_time_1d < end_time_1d  .AND.  .NOT. stop_dt_1d )
380
381!
382!--    Depending on the timestep scheme, carry out one or more intermediate
383!--    timesteps
384
385       intermediate_timestep_count = 0
386       DO  WHILE ( intermediate_timestep_count < &
387                   intermediate_timestep_count_max )
388
389          intermediate_timestep_count = intermediate_timestep_count + 1
390
391          CALL timestep_scheme_steering
392
393!
394!--       Compute all tendency terms. If a Prandtl-layer is simulated, k starts
395!--       at nzb+2.
396          DO  k = nzb_diff, nzt
397
398             kmzm = 0.5_wp * ( km1d(k-1) + km1d(k) )
399             kmzp = 0.5_wp * ( km1d(k) + km1d(k+1) )
400!
401!--          u-component
402             te_u(k) =  f * ( v1d(k) - vg(k) ) + ( &
403                              kmzp * ( u1d(k+1) - u1d(k) ) * ddzu(k+1) &
404                            - kmzm * ( u1d(k) - u1d(k-1) ) * ddzu(k)   &
405                                                 ) * ddzw(k)
406!
407!--          v-component
408             te_v(k) = -f * ( u1d(k) - ug(k) ) + (                     &
409                              kmzp * ( v1d(k+1) - v1d(k) ) * ddzu(k+1) &
410                            - kmzm * ( v1d(k) - v1d(k-1) ) * ddzu(k)   &
411                                                 ) * ddzw(k)
412          ENDDO
413          IF ( .NOT. constant_diffusion )  THEN
414             DO  k = nzb_diff, nzt
415!
416!--             TKE
417                kmzm = 0.5_wp * ( km1d(k-1) + km1d(k) )
418                kmzp = 0.5_wp * ( km1d(k) + km1d(k+1) )
419                IF ( .NOT. humidity )  THEN
420                   pt_0 = pt_init(k)
421                   flux =  ( pt_init(k+1)-pt_init(k-1) ) * dd2zu(k)
422                ELSE
423                   pt_0 = pt_init(k) * ( 1.0_wp + 0.61_wp * q_init(k) )
424                   flux = ( ( pt_init(k+1) - pt_init(k-1) ) +                  &
425                            0.61_wp * ( pt_init(k+1) * q_init(k+1) -           &
426                                        pt_init(k-1) * q_init(k-1)   )         &
427                          ) * dd2zu(k)
428                ENDIF
429
430                IF ( dissipation_1d == 'detering' )  THEN
431!
432!--                According to Detering, c_e=c_m^3
433                   dissipation = c_m**3 * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
434                ELSEIF ( dissipation_1d == 'as_in_3d_model' )  THEN
435                   dissipation = ( 0.19_wp                                     &
436                                   + 0.74_wp * l1d_diss(k) /  l_grid(k)        &
437                                 ) * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
438                ENDIF
439
440                te_e(k) = km1d(k) * ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2&
441                                    + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2&
442                                    )                                          &
443                                    - g / pt_0 * kh1d(k) * flux                &
444                                    +            (                             &
445                                     kmzp * ( e1d(k+1) - e1d(k) ) * ddzu(k+1)  &
446                                   - kmzm * ( e1d(k) - e1d(k-1) ) * ddzu(k)    &
447                                                 ) * ddzw(k)                   &
448                                   - dissipation
449             ENDDO
450          ENDIF
451
452!
453!--       Tendency terms at the top of the Prandtl-layer.
454!--       Finite differences of the momentum fluxes are computed using half the
455!--       normal grid length (2.0*ddzw(k)) for the sake of enhanced accuracy
456          IF ( constant_flux_layer )  THEN
457
458             k = nzb+1
459             kmzm = 0.5_wp * ( km1d(k-1) + km1d(k) )
460             kmzp = 0.5_wp * ( km1d(k) + km1d(k+1) )
461             IF ( .NOT. humidity )  THEN
462                pt_0 = pt_init(k)
463                flux =  ( pt_init(k+1)-pt_init(k-1) ) * dd2zu(k)
464             ELSE
465                pt_0 = pt_init(k) * ( 1.0_wp + 0.61_wp * q_init(k) )
466                flux = ( ( pt_init(k+1) - pt_init(k-1) ) +                     &
467                         0.61_wp * ( pt_init(k+1) * q_init(k+1) -              &
468                                     pt_init(k-1) * q_init(k-1)   )            &
469                       ) * dd2zu(k)
470             ENDIF
471
472             IF ( dissipation_1d == 'detering' )  THEN
473!
474!--             According to Detering, c_e=0.064
475                dissipation = 0.064_wp * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
476             ELSEIF ( dissipation_1d == 'as_in_3d_model' )  THEN
477                dissipation = ( 0.19_wp + 0.74_wp * l1d_diss(k) / l_grid(k) )  &
478                              * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
479             ENDIF
480
481!
482!--          u-component
483             te_u(k) = f * ( v1d(k) - vg(k) ) + (                              &
484                       kmzp * ( u1d(k+1) - u1d(k) ) * ddzu(k+1) + usws1d       &
485                                                ) * 2.0_wp * ddzw(k)
486!
487!--          v-component
488             te_v(k) = -f * ( u1d(k) - ug(k) ) + (                             &
489                       kmzp * ( v1d(k+1) - v1d(k) ) * ddzu(k+1) + vsws1d       &
490                                                 ) * 2.0_wp * ddzw(k)
491!
492!--          TKE
493             te_e(k) = km1d(k) * ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2   &
494                                 + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2   &
495                                 )                                             &
496                                 - g / pt_0 * kh1d(k) * flux                   &
497                                 +           (                                 &
498                                  kmzp * ( e1d(k+1) - e1d(k) ) * ddzu(k+1)     &
499                                - kmzm * ( e1d(k) - e1d(k-1) ) * ddzu(k)       &
500                                              ) * ddzw(k)                      &
501                                - dissipation
502          ENDIF
503
504!
505!--       Prognostic equations for all 1D variables
506          DO  k = nzb+1, nzt
507
508             u1d_p(k) = u1d(k) + dt_1d * ( tsc(2) * te_u(k) + &
509                                           tsc(3) * te_um(k) )
510             v1d_p(k) = v1d(k) + dt_1d * ( tsc(2) * te_v(k) + &
511                                           tsc(3) * te_vm(k) )
512
513          ENDDO
514          IF ( .NOT. constant_diffusion )  THEN
515             DO  k = nzb+1, nzt
516
517                e1d_p(k) = e1d(k) + dt_1d * ( tsc(2) * te_e(k) + &
518                                              tsc(3) * te_em(k) )
519
520             ENDDO
521!
522!--          Eliminate negative TKE values, which can result from the
523!--          integration due to numerical inaccuracies. In such cases the TKE
524!--          value is reduced to 10 percent of its old value.
525             WHERE ( e1d_p < 0.0_wp )  e1d_p = 0.1_wp * e1d
526          ENDIF
527
528!
529!--       Calculate tendencies for the next Runge-Kutta step
530          IF ( timestep_scheme(1:5) == 'runge' ) THEN
531             IF ( intermediate_timestep_count == 1 )  THEN
532
533                DO  k = nzb+1, nzt
534                   te_um(k) = te_u(k)
535                   te_vm(k) = te_v(k)
536                ENDDO
537
538                IF ( .NOT. constant_diffusion )  THEN
539                   DO k = nzb+1, nzt
540                      te_em(k) = te_e(k)
541                   ENDDO
542                ENDIF
543
544             ELSEIF ( intermediate_timestep_count < &
545                         intermediate_timestep_count_max )  THEN
546
547                DO  k = nzb+1, nzt
548                   te_um(k) = -9.5625_wp * te_u(k) + 5.3125_wp * te_um(k)
549                   te_vm(k) = -9.5625_wp * te_v(k) + 5.3125_wp * te_vm(k)
550                ENDDO
551
552                IF ( .NOT. constant_diffusion )  THEN
553                   DO k = nzb+1, nzt
554                      te_em(k) = -9.5625_wp * te_e(k) + 5.3125_wp * te_em(k)
555                   ENDDO
556                ENDIF
557
558             ENDIF
559          ENDIF
560
561
562!
563!--       Boundary conditions for the prognostic variables.
564!--       At the top boundary (nzt+1) u,v and e keep their initial values
565!--       (ug(nzt+1), vg(nzt+1), 0).
566!--       At the bottom boundary, Dirichlet condition is used for u and v (0)
567!--       and Neumann condition for e (e(nzb)=e(nzb+1)).
568          u1d_p(nzb) = 0.0_wp
569          v1d_p(nzb) = 0.0_wp
570
571!
572!--       Swap the time levels in preparation for the next time step.
573          u1d  = u1d_p
574          v1d  = v1d_p
575          IF ( .NOT. constant_diffusion )  THEN
576             e1d  = e1d_p
577          ENDIF
578
579!
580!--       Compute diffusion quantities
581          IF ( .NOT. constant_diffusion )  THEN
582
583!
584!--          First compute the vertical fluxes in the Prandtl-layer
585             IF ( constant_flux_layer )  THEN
586!
587!--             Compute theta* using Rif numbers of the previous time step
588                IF ( rif1d(nzb+1) >= 0.0_wp )  THEN
589!
590!--                Stable stratification
591                   ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) /          &
592                          ( LOG( zu(nzb+1) / z0h1d ) + 5.0_wp * rif1d(nzb+1) * &
593                                          ( zu(nzb+1) - z0h1d ) / zu(nzb+1)    &
594                          )
595                ELSE
596!
597!--                Unstable stratification
598                   a = SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) )
599                   b = SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) /                 &
600                       zu(nzb+1) * z0h1d )
601
602                   ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) /          &
603                          LOG( (a-1.0_wp) / (a+1.0_wp) *                       &
604                               (b+1.0_wp) / (b-1.0_wp) )
605                ENDIF
606
607             ENDIF    ! constant_flux_layer
608
609!
610!--          Compute the Richardson-flux numbers,
611!--          first at the top of the Prandtl-layer using u* of the previous
612!--          time step (+1E-30, if u* = 0), then in the remaining area. There
613!--          the rif-numbers of the previous time step are used.
614
615             IF ( constant_flux_layer )  THEN
616                IF ( .NOT. humidity )  THEN
617                   pt_0 = pt_init(nzb+1)
618                   flux = ts1d
619                ELSE
620                   pt_0 = pt_init(nzb+1) * ( 1.0_wp + 0.61_wp * q_init(nzb+1) )
621                   flux = ts1d + 0.61_wp * pt_init(k) * qs1d
622                ENDIF
623                rif1d(nzb+1) = zu(nzb+1) * kappa * g * flux / &
624                               ( pt_0 * ( us1d**2 + 1E-30_wp ) )
625             ENDIF
626
627             DO  k = nzb_diff, nzt
628                IF ( .NOT. humidity )  THEN
629                   pt_0 = pt_init(k)
630                   flux = ( pt_init(k+1) - pt_init(k-1) ) * dd2zu(k)
631                ELSE
632                   pt_0 = pt_init(k) * ( 1.0_wp + 0.61_wp * q_init(k) )
633                   flux = ( ( pt_init(k+1) - pt_init(k-1) )                    &
634                            + 0.61_wp                                          &
635                            * (   pt_init(k+1) * q_init(k+1)                   &
636                                - pt_init(k-1) * q_init(k-1) )                 &
637                          ) * dd2zu(k)
638                ENDIF
639                IF ( rif1d(k) >= 0.0_wp )  THEN
640                   rif1d(k) = g / pt_0 * flux /                                &
641                              (  ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2     &
642                               + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2     &
643                               + 1E-30_wp                                      &
644                              )
645                ELSE
646                   rif1d(k) = g / pt_0 * flux /                                &
647                              (  ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2     &
648                               + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2     &
649                               + 1E-30_wp                                      &
650                              ) * ( 1.0_wp - 16.0_wp * rif1d(k) )**0.25_wp
651                ENDIF
652             ENDDO
653!
654!--          Richardson-numbers must remain restricted to a realistic value
655!--          range. It is exceeded excessively for very small velocities
656!--          (u,v --> 0).
657             WHERE ( rif1d < -5.0_wp )  rif1d = -5.0_wp
658             WHERE ( rif1d > 1.0_wp )  rif1d = 1.0_wp
659
660!
661!--          Compute u* from the absolute velocity value
662             IF ( constant_flux_layer )  THEN
663                uv_total = SQRT( u1d(nzb+1)**2 + v1d(nzb+1)**2 )
664
665                IF ( rif1d(nzb+1) >= 0.0_wp )  THEN
666!
667!--                Stable stratification
668                   us1d = kappa * uv_total / (                                 &
669                             LOG( zu(nzb+1) / z01d ) + 5.0_wp * rif1d(nzb+1) * &
670                                              ( zu(nzb+1) - z01d ) / zu(nzb+1) &
671                                             )
672                ELSE
673!
674!--                Unstable stratification
675                   a = 1.0_wp / SQRT( SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) ) )
676                   b = 1.0_wp / SQRT( SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) /  &
677                                                     zu(nzb+1) * z01d ) )
678                   us1d = kappa * uv_total / (                                 &
679                              LOG( (1.0_wp+b) / (1.0_wp-b) * (1.0_wp-a) /      &
680                                   (1.0_wp+a) ) +                              &
681                              2.0_wp * ( ATAN( b ) - ATAN( a ) )               &
682                                             )
683                ENDIF
684
685!
686!--             Compute the momentum fluxes for the diffusion terms
687                usws1d  = - u1d(nzb+1) / uv_total * us1d**2
688                vsws1d  = - v1d(nzb+1) / uv_total * us1d**2
689
690!
691!--             Boundary condition for the turbulent kinetic energy at the top
692!--             of the Prandtl-layer. c_m = 0.4 according to Detering.
693!--             Additional Neumann condition de/dz = 0 at nzb is set to ensure
694!--             compatibility with the 3D model.
695                IF ( ibc_e_b == 2 )  THEN
696                   e1d(nzb+1) = ( us1d / c_m )**2
697                ENDIF
698                e1d(nzb) = e1d(nzb+1)
699
700                IF ( humidity ) THEN
701!
702!--                Compute q*
703                   IF ( rif1d(nzb+1) >= 0.0_wp )  THEN
704!
705!--                   Stable stratification
706                      qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) /         &
707                          ( LOG( zu(nzb+1) / z0h1d ) + 5.0_wp * rif1d(nzb+1) * &
708                                          ( zu(nzb+1) - z0h1d ) / zu(nzb+1)    &
709                          )
710                   ELSE
711!
712!--                   Unstable stratification
713                      a = SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) )
714                      b = SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) /              &
715                                         zu(nzb+1) * z0h1d )
716                      qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) /         &
717                             LOG( (a-1.0_wp) / (a+1.0_wp) *                    &
718                                  (b+1.0_wp) / (b-1.0_wp) )
719                   ENDIF
720                ELSE
721                   qs1d = 0.0_wp
722                ENDIF
723
724             ENDIF   !  constant_flux_layer
725
726!
727!--          Compute the diabatic mixing length. The unstable stratification
728!--          must not be considered for l1d (km1d) as it is already considered
729!--          in the dissipation of TKE via l1d_diss. Otherwise, km1d would be
730!--          too large.
731             IF ( mixing_length_1d == 'blackadar' )  THEN
732                DO  k = nzb+1, nzt
733                   IF ( rif1d(k) >= 0.0_wp )  THEN
734                      l1d(k) = l_black(k) / ( 1.0_wp + 5.0_wp * rif1d(k) )
735                      l1d_diss(k) = l1d(k)
736                   ELSE
737                      l1d(k) = l_black(k)
738                      l1d_diss(k) = l_black(k) *                               &
739                                    SQRT( 1.0_wp - 16.0_wp * rif1d(k) )
740                   ENDIF
741                ENDDO
742             ELSEIF ( mixing_length_1d == 'as_in_3d_model' )  THEN
743                DO  k = nzb+1, nzt
744                   dpt_dz = ( pt_init(k+1) - pt_init(k-1) ) * dd2zu(k)
745                   IF ( dpt_dz > 0.0_wp )  THEN
746                      l_stable = 0.76_wp * SQRT( e1d(k) ) /                    &
747                                     SQRT( g / pt_init(k) * dpt_dz ) + 1E-5_wp
748                   ELSE
749                      l_stable = l_grid(k)
750                   ENDIF
751                   l1d(k) = MIN( l_grid(k), l_stable )
752                   l1d_diss(k) = l1d(k)
753                ENDDO
754             ENDIF
755
756!
757!--          Compute the diffusion coefficients for momentum via the
758!--          corresponding Prandtl-layer relationship and according to
759!--          Prandtl-Kolmogorov, respectively
760             IF ( constant_flux_layer )  THEN
761                IF ( rif1d(nzb+1) >= 0.0_wp )  THEN
762                   km1d(nzb+1) = us1d * kappa * zu(nzb+1) /                    &
763                                 ( 1.0_wp + 5.0_wp * rif1d(nzb+1) )
764                ELSE
765                   km1d(nzb+1) = us1d * kappa * zu(nzb+1) *                    &
766                                 ( 1.0_wp - 16.0_wp * rif1d(nzb+1) )**0.25_wp
767                ENDIF
768             ENDIF
769             DO  k = nzb_diff, nzt
770                km1d(k) = c_m * SQRT( e1d(k) ) * l1d(k)
771             ENDDO
772
773!
774!--          Add damping layer
775             DO  k = damp_level_ind_1d+1, nzt+1
776                km1d(k) = 1.1_wp * km1d(k-1)
777                km1d(k) = MIN( km1d(k), 10.0_wp )
778             ENDDO
779
780!
781!--          Compute the diffusion coefficient for heat via the relationship
782!--          kh = phim / phih * km
783             DO  k = nzb+1, nzt
784                IF ( rif1d(k) >= 0.0_wp )  THEN
785                   kh1d(k) = km1d(k)
786                ELSE
787                   kh1d(k) = km1d(k) * ( 1.0_wp - 16.0_wp * rif1d(k) )**0.25_wp
788                ENDIF
789             ENDDO
790
791          ENDIF   ! .NOT. constant_diffusion
792
793       ENDDO   ! intermediate step loop
794
795!
796!--    Increment simulated time and output times
797       current_timestep_number_1d = current_timestep_number_1d + 1
798       simulated_time_1d          = simulated_time_1d + dt_1d
799       simulated_time_chr         = time_to_string( simulated_time_1d )
800       time_pr_1d                 = time_pr_1d          + dt_1d
801       time_run_control_1d        = time_run_control_1d + dt_1d
802
803!
804!--    Determine and print out quantities for run control
805       IF ( time_run_control_1d >= dt_run_control_1d )  THEN
806          CALL run_control_1d
807          time_run_control_1d = time_run_control_1d - dt_run_control_1d
808       ENDIF
809
810!
811!--    Profile output on file
812       IF ( time_pr_1d >= dt_pr_1d )  THEN
813          CALL print_1d_model
814          time_pr_1d = time_pr_1d - dt_pr_1d
815       ENDIF
816
817!
818!--    Determine size of next time step
819       CALL timestep_1d
820
821    ENDDO   ! time loop
822
823
824 END SUBROUTINE time_integration_1d
825
826
827!------------------------------------------------------------------------------!
828! Description:
829! ------------
830!> Compute and print out quantities for run control of the 1D model.
831!------------------------------------------------------------------------------!
832 
833 SUBROUTINE run_control_1d
834
835
836    USE constants,                                                             &
837        ONLY:  pi
838
839    IMPLICIT NONE
840
841    INTEGER(iwp) ::  k     !< loop index
842   
843    REAL(wp) ::  alpha     !< angle of wind vector at top of constant-flux layer
844    REAL(wp) ::  energy    !< kinetic energy
845    REAL(wp) ::  umax      !< maximum of u
846    REAL(wp) ::  uv_total  !< horizontal wind speed
847    REAL(wp) ::  vmax      !< maximum of v
848
849!
850!-- Output
851    IF ( myid == 0 )  THEN
852!
853!--    If necessary, write header
854       IF ( .NOT. run_control_header_1d )  THEN
855          CALL check_open( 15 )
856          WRITE ( 15, 100 )
857          run_control_header_1d = .TRUE.
858       ENDIF
859
860!
861!--    Compute control quantities
862!--    grid level nzp is excluded due to mirror boundary condition
863       umax = 0.0_wp; vmax = 0.0_wp; energy = 0.0_wp
864       DO  k = nzb+1, nzt+1
865          umax = MAX( ABS( umax ), ABS( u1d(k) ) )
866          vmax = MAX( ABS( vmax ), ABS( v1d(k) ) )
867          energy = energy + 0.5_wp * ( u1d(k)**2 + v1d(k)**2 )
868       ENDDO
869       energy = energy / REAL( nzt - nzb + 1, KIND=wp )
870
871       uv_total = SQRT( u1d(nzb+1)**2 + v1d(nzb+1)**2 )
872       IF ( ABS( v1d(nzb+1) ) < 1.0E-5_wp )  THEN
873          alpha = ACOS( SIGN( 1.0_wp , u1d(nzb+1) ) )
874       ELSE
875          alpha = ACOS( u1d(nzb+1) / uv_total )
876          IF ( v1d(nzb+1) <= 0.0_wp )  alpha = 2.0_wp * pi - alpha
877       ENDIF
878       alpha = alpha / ( 2.0_wp * pi ) * 360.0_wp
879
880       WRITE ( 15, 101 )  current_timestep_number_1d, simulated_time_chr, &
881                          dt_1d, umax, vmax, us1d, alpha, energy
882!
883!--    Write buffer contents to disc immediately
884       FLUSH( 15 )
885
886    ENDIF
887
888!
889!-- formats
890100 FORMAT (///'1D run control output:'/ &
891              &'------------------------------'// &
892           &'ITER.  HH:MM:SS    DT      UMAX   VMAX    U*   ALPHA   ENERG.'/ &
893           &'-------------------------------------------------------------')
894101 FORMAT (I5,2X,A9,1X,F6.2,2X,F6.2,1X,F6.2,1X,F6.3,2X,F5.1,2X,F7.2)
895
896
897 END SUBROUTINE run_control_1d
898
899
900
901!------------------------------------------------------------------------------!
902! Description:
903! ------------
904!> Compute the time step w.r.t. the diffusion criterion
905!------------------------------------------------------------------------------!
906 
907 SUBROUTINE timestep_1d
908
909    IMPLICIT NONE
910
911    INTEGER(iwp) ::  k    !< loop index
912
913    REAL(wp) ::  dt_diff  !< time step accorind to diffusion criterion
914    REAL(wp) ::  fac      !< factor of criterion
915    REAL(wp) ::  value    !< auxiliary variable
916
917!
918!-- Compute the currently feasible time step according to the diffusion
919!-- criterion. At nzb+1 the half grid length is used.
920    fac = 0.35_wp
921    dt_diff = dt_max_1d
922    DO  k = nzb+2, nzt
923       value   = fac * dzu(k) * dzu(k) / ( km1d(k) + 1E-20_wp )
924       dt_diff = MIN( value, dt_diff )
925    ENDDO
926    value   = fac * zu(nzb+1) * zu(nzb+1) / ( km1d(nzb+1) + 1E-20_wp )
927    dt_1d = MIN( value, dt_diff )
928
929!
930!-- Set flag when the time step becomes too small
931    IF ( dt_1d < ( 0.00001_wp * dt_max_1d ) )  THEN
932       stop_dt_1d = .TRUE.
933
934       WRITE( message_string, * ) 'timestep has exceeded the lower limit &', &
935                                  'dt_1d = ',dt_1d,' s   simulation stopped!'
936       CALL message( 'timestep_1d', 'PA0192', 1, 2, 0, 6, 0 )
937       
938    ENDIF
939
940 END SUBROUTINE timestep_1d
941
942
943
944!------------------------------------------------------------------------------!
945! Description:
946! ------------
947!> List output of profiles from the 1D-model
948!------------------------------------------------------------------------------!
949 
950 SUBROUTINE print_1d_model
951
952    IMPLICIT NONE
953
954    INTEGER(iwp) ::  k  !< loop parameter
955
956    LOGICAL, SAVE :: write_first = .TRUE. !< flag for writing header
957
958
959    IF ( myid == 0 )  THEN
960!
961!--    Open list output file for profiles from the 1D-model
962       CALL check_open( 17 )
963
964!
965!--    Write Header
966       IF ( write_first )  THEN
967          WRITE ( 17, 100 )  TRIM( run_description_header )
968          write_first = .FALSE.
969       ENDIF
970
971!
972!--    Write the values
973       WRITE ( 17, 104 )  TRIM( simulated_time_chr )
974       WRITE ( 17, 101 )
975       WRITE ( 17, 102 )
976       WRITE ( 17, 101 )
977       DO  k = nzt+1, nzb, -1
978          WRITE ( 17, 103)  k, zu(k), u1d(k), v1d(k), pt_init(k), e1d(k), &
979                            rif1d(k), km1d(k), kh1d(k), l1d(k), zu(k), k
980       ENDDO
981       WRITE ( 17, 101 )
982       WRITE ( 17, 102 )
983       WRITE ( 17, 101 )
984
985!
986!--    Write buffer contents to disc immediately
987       FLUSH( 17 )
988
989    ENDIF
990
991!
992!-- Formats
993100 FORMAT ('# ',A/'#',10('-')/'# 1d-model profiles')
994104 FORMAT (//'# Time: ',A)
995101 FORMAT ('#',79('-'))
996102 FORMAT ('#  k     zu      u      v     pt      e    rif    Km    Kh     ', &
997            'l      zu      k')
998103 FORMAT (1X,I4,1X,F7.1,1X,F6.2,1X,F6.2,1X,F6.2,1X,F6.2,1X,F5.2,1X,F5.2, &
999            1X,F5.2,1X,F6.2,1X,F7.1,2X,I4)
1000
1001
1002 END SUBROUTINE print_1d_model
1003
1004
1005 END MODULE
Note: See TracBrowser for help on using the repository browser.