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

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