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

Last change on this file since 1 was 1, checked in by raasch, 15 years ago

Initial repository layout and content

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