source: palm/trunk/SOURCE/init_1d_model.f90 @ 2320

Last change on this file since 2320 was 2299, checked in by maronga, 7 years ago

improvements for spinup mechanism

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