source: palm/tags/release-3.2/SOURCE/init_1d_model.f90 @ 4109

Last change on this file since 4109 was 77, checked in by raasch, 17 years ago

New:
---

particle reflection from vertical walls implemented, particle SGS model adjusted to walls

Wall functions for vertical walls now include diabatic conditions. New subroutines wall_fluxes, wall_fluxes_e. New 4D-array rif_wall.

new d3par-parameter netcdf_64bit_3d to switch on 64bit offset only for 3D files

new d3par-parameter dt_max to define the maximum value for the allowed timestep

new inipar-parameter loop_optimization to control the loop optimization method

new inipar-parameter pt_refrence. If given, this value is used as the reference that in buoyancy terms (otherwise, the instantaneous horizontally averaged temperature is used).

new user interface user_advec_particles

new initializing action "by_user" calls user_init_3d_model and allows the initial setting of all 3d arrays

topography height informations are stored on arrays zu_s_inner and zw_w_inner and output to the 2d/3d NetCDF files

samples added to the user interface which show how to add user-define time series quantities.

calculation/output of precipitation amount, precipitation rate and z0 (by setting "pra*", "prr*", "z0*" with data_output). The time interval on which the precipitation amount is defined is set by new d3par-parameter precipitation_amount_interval

unit 9 opened for debug output (file DEBUG_<pe#>)

Makefile, advec_particles, average_3d_data, buoyancy, calc_precipitation, check_open, check_parameters, data_output_2d, diffusion_e, diffusion_u, diffusion_v, diffusion_w, diffusivities, header, impact_of_latent_heat, init_particles, init_3d_model, modules, netcdf, parin, production_e, read_var_list, read_3d_binary, sum_up_3d_data, user_interface, write_var_list, write_3d_binary

New: wall_fluxes

Changed:


General revision of non-cyclic horizontal boundary conditions: radiation boundary conditions are now used instead of Neumann conditions at the outflow (calculation needs velocity values for t-dt, which are stored on new arrays u_m_l, u_m_r, etc.), calculation of mean outflow is not needed any more, volume flow control is added for the outflow boundary (currently only for the north boundary!!), additional gridpoints along x and y (uxrp, vynp) are not needed any more, routine "boundary_conds" now operates on timelevel t+dt and is not split in two parts (main, uvw_outflow) any more, Neumann boundary conditions at inflow/outflow in case of non-cyclic boundary conditions for all 2d-arrays that are handled by exchange_horiz_2d

The FFT-method for solving the Poisson-equation is now working with Neumann boundary conditions both at the bottom and the top. This requires adjustments of the tridiagonal coefficients and subtracting the horizontally averaged mean from the vertical velocity field.

+age_m in particle_type

Particles-package is now part of the default code ("-p particles" is not needed any more).

Move call of user_actions( 'after_integration' ) below increment of times
and counters. user_actions is now called for each statistic region and has as an argument the number of the respective region (sr)

d3par-parameter data_output_ts removed. Timeseries output for "profil" removed. Timeseries are now switched on by dt_dots. Timeseries data is collected in flow_statistics.

Initial velocities at nzb+1 are regarded for volume flow control in case they have been set zero before (to avoid small timesteps); see new internal parameters u/v_nzb_p1_for_vfc.

q is not allowed to become negative (prognostic_equations).

poisfft_init is only called if fft-solver is switched on (init_pegrid).

d3par-parameter moisture renamed to humidity.

Subversion global revision number is read from mrun and added to the run description header and to the run control (_rc) file.

vtk directives removed from main program.

The uitility routine interpret_config reads PALM environment variables from NAMELIST instead using the system call GETENV.

advec_u_pw, advec_u_up, advec_v_pw, advec_v_up, asselin_filter, check_parameters, coriolis, data_output_dvrp, data_output_ptseries, data_output_ts, data_output_2d, data_output_3d, diffusion_u, diffusion_v, exchange_horiz, exchange_horiz_2d, flow_statistics, header, init_grid, init_particles, init_pegrid, init_rankine, init_pt_anomaly, init_1d_model, init_3d_model, modules, palm, package_parin, parin, poisfft, poismg, prandtl_fluxes, pres, production_e, prognostic_equations, read_var_list, read_3d_binary, sor, swap_timelevel, time_integration, write_var_list, write_3d_binary

Errors:


Bugfix: preset of tendencies te_em, te_um, te_vm in init_1d_model

Bugfix in sample for reading user defined data from restart file (user_init)

Bugfix in setting diffusivities for cases with the outflow damping layer extending over more than one subdomain (init_3d_model)

Check for possible negative humidities in the initial humidity profile.

in Makefile, default suffixes removed from the suffix list to avoid calling of m2c in
# case of .mod files

Makefile
check_parameters, init_1d_model, init_3d_model, user_interface

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