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

Last change on this file since 76 was 75, checked in by raasch, 17 years ago

preliminary update for changes concerning non-cyclic boundary conditions

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