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

Last change on this file since 1322 was 1322, checked in by raasch, 8 years ago

REAL functions and a lot of REAL constants provided with KIND-attribute,
some small bugfixes

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