source: palm/tags/release-3.1b/SOURCE/init_1d_model.f90 @ 226

Last change on this file since 226 was 4, checked in by raasch, 17 years ago

Id keyword set as property for all *.f90 files

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