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

Last change on this file since 1320 was 1320, checked in by raasch, 7 years ago

ONLY-attribute added to USE-statements,
kind-parameters added to all INTEGER and REAL declaration statements,
kinds are defined in new module kinds,
old module precision_kind is removed,
revision history before 2012 removed,
comment fields (!:) to be used for variable explanations added to all variable declaration statements

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