source: palm/tags/release-3.4/SOURCE/init_1d_model.f90 @ 3911

Last change on this file since 3911 was 83, checked in by raasch, 17 years ago

New:
---

Changed:


PALM can be generally installed on any kind of Linux-, IBM-AIX-, or NEC-SX-system by adding appropriate settings to the configuration file.

Scripts are also running under the public domain ksh.

All system relevant compile and link options as well as the host identifier (local_host) are specified in the configuration file.

Filetransfer by ftp removed (options -f removed from mrun and mbuild).

Call of (system-)FLUSH routine moved to new routine local_flush.

return_addres and return_username are read from ENVPAR-NAMELIST-file instead of using local_getenv.

Preprocessor strings for different linux clusters changed to "lc", some preprocessor directives renamed (new: intel_openmp_bug), preprocessor directives for old systems removed

advec_particles, check_open, cpu_log, cpu_statistics, data_output_dvrp, flow_statistics, header, init_dvrp, init_particles, init_1d_model, init_dvrp, init_pegrid, local_getenv, local_system, local_tremain, local_tremain_ini, modules, palm, parin, run_control

new:
local_flush

mbuild, mrun

Errors:


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