source: palm/trunk/SOURCE/check_parameters.f90 @ 87

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

Preliminary update for user defined profiles

  • Property svn:keywords set to Id
File size: 89.4 KB
RevLine 
[1]1 SUBROUTINE check_parameters
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
[87]6! Check for user-defined profiles
[77]7!
8! Former revisions:
9! -----------------
10! $Id: check_parameters.f90 87 2007-05-22 15:46:47Z raasch $
11!
12! 75 2007-03-22 09:54:05Z raasch
[51]13! "by_user" allowed as initializing action, -data_output_ts,
[63]14! leapfrog with non-flat topography not allowed any more, loop_optimization
[75]15! and pt_reference are checked, moisture renamed humidity,
[72]16! output of precipitation amount/rate and roughnes length + check
[73]17! possible negative humidities are avoided in initial profile,
[75]18! dirichlet/neumann changed to dirichlet/radiation, etc.,
19! revision added to run_description_header
[1]20!
[39]21! 20 2007-02-26 00:12:32Z raasch
22! Temperature and humidity gradients at top are now calculated for nzt+1,
23! top_heatflux and respective boundary condition bc_pt_t is checked
24!
[3]25! RCS Log replace by Id keyword, revision history cleaned up
26!
[1]27! Revision 1.61  2006/08/04 14:20:25  raasch
28! do2d_unit and do3d_unit now defined as 2d-arrays, check of
29! use_upstream_for_tke, default value for dt_dopts,
30! generation of file header moved from routines palm and header to here
31!
32! Revision 1.1  1997/08/26 06:29:23  raasch
33! Initial revision
34!
35!
36! Description:
37! ------------
38! Check control parameters and deduce further quantities.
39!------------------------------------------------------------------------------!
40
41    USE arrays_3d
42    USE constants
43    USE control_parameters
44    USE grid_variables
45    USE indices
46    USE model_1d
47    USE netcdf_control
48    USE particle_attributes
49    USE pegrid
50    USE profil_parameter
51    USE statistics
52    USE transpose_indices
53
54    IMPLICIT NONE
55
56    CHARACTER (LEN=1)   ::  sq
57    CHARACTER (LEN=6)   ::  var
58    CHARACTER (LEN=7)   ::  unit
59    CHARACTER (LEN=8)   ::  date
60    CHARACTER (LEN=10)  ::  time
61    CHARACTER (LEN=100) ::  action
62
63    INTEGER ::  i, ilen, intervals, iter, j, k, nnxh, nnyh, position, prec
64    LOGICAL ::  found, ldum
65    REAL    ::  gradient, maxn, maxp
66
67
68!
69!-- Warning, if host is not set
70    IF ( host(1:1) == ' ' )  THEN
71       IF ( myid == 0 )  THEN
72          PRINT*, '+++ WARNING: check_parameters:'
73          PRINT*, '    "host" is not set. Please check that environment', &
74                       ' variable "localhost"'
75          PRINT*, '    is set before running PALM'
76       ENDIF
77    ENDIF
78
79!
80!-- Generate the file header which is used as a header for most of PALM's
81!-- output files
82    CALL DATE_AND_TIME( date, time )
83    run_date = date(7:8)//'-'//date(5:6)//'-'//date(3:4)
84    run_time = time(1:2)//':'//time(3:4)//':'//time(5:6)
85
[75]86    WRITE ( run_description_header, '(A,2X,A,2X,A,A,A,I2.2,2X,A,A,2X,A,1X,A)' )&
87              TRIM( version ), TRIM( revision ), &
88              'run: ', TRIM( run_identifier ), '.', &
[1]89              runnr, 'host: ', TRIM( host ), run_date, run_time
90
91!
[63]92!-- Check the general loop optimization method
93    IF ( loop_optimization == 'default' )  THEN
94       IF ( host(1:3) == 'nec' )  THEN
95          loop_optimization = 'vector'
96       ELSE
97          loop_optimization = 'cache'
98       ENDIF
99    ENDIF
100    IF ( loop_optimization /= 'noopt'  .AND.  loop_optimization /= 'cache' &
101         .AND.  loop_optimization /= 'vector' )  THEN
102       IF ( myid == 0 )  THEN
103          PRINT*, '+++ check_parameters:'
104          PRINT*, '    illegal value given for loop_optimization: ', &
105                  TRIM( loop_optimization )
106       ENDIF
107       CALL local_stop
108    ENDIF
109
110!
[1]111!-- Check topography setting (check for illegal parameter combinations)
112    IF ( topography /= 'flat' )  THEN
113       action = ' '
114       IF ( scalar_advec /= 'pw-scheme' )  THEN
115          WRITE( action, '(A,A)' )  'scalar_advec = ', scalar_advec
116       ENDIF
117       IF ( momentum_advec /= 'pw-scheme' )  THEN
118          WRITE( action, '(A,A)' )  'momentum_advec = ', momentum_advec
119       ENDIF
[51]120       IF ( timestep_scheme(1:8) == 'leapfrog' )  THEN
121          WRITE( action, '(A,A)' )  'timestep_scheme = ', timestep_scheme
122       ENDIF
[1]123       IF ( psolver == 'multigrid'  .OR.  psolver == 'sor' )  THEN
124          WRITE( action, '(A,A)' )  'psolver = ', psolver
125       ENDIF
126       IF ( sloping_surface )  THEN
127          WRITE( action, '(A)' )  'sloping surface = .TRUE.'
128       ENDIF
129       IF ( galilei_transformation )  THEN
130          WRITE( action, '(A)' )  'galilei_transformation = .TRUE.'
131       ENDIF
132       IF ( cloud_physics )  THEN
133          WRITE( action, '(A)' )  'cloud_physics = .TRUE.'
134       ENDIF
135       IF ( cloud_droplets )  THEN
136          WRITE( action, '(A)' )  'cloud_droplets = .TRUE.'
137       ENDIF
[75]138       IF ( humidity )  THEN
139          WRITE( action, '(A)' )  'humidity = .TRUE.'
[1]140       ENDIF
141       IF ( .NOT. prandtl_layer )  THEN
142          WRITE( action, '(A)' )  'prandtl_layer = .FALSE.'
143       ENDIF
144       IF ( action /= ' ' )  THEN
145          IF ( myid == 0 )  THEN
146             PRINT*, '+++ check_parameters:'
147             PRINT*, '    a non-flat topography does not allow ', TRIM( action )
148          ENDIF
149          CALL local_stop
150       ENDIF
151    ENDIF
152!
153!-- Check whether there are any illegal values
154!-- Pressure solver:
155    IF ( psolver /= 'poisfft'  .AND.  psolver /= 'poisfft_hybrid'  .AND. &
156         psolver /= 'sor'  .AND.  psolver /= 'multigrid' )  THEN
157       IF ( myid == 0 )  THEN
158          PRINT*, '+++ check_parameters:'
159          PRINT*, '    unknown solver for perturbation pressure: psolver=', &
160                  psolver
161       ENDIF
162       CALL local_stop
163    ENDIF
164
165#if defined( __parallel )
166    IF ( psolver == 'poisfft_hybrid'  .AND.  pdims(2) /= 1 )  THEN
167       IF ( myid == 0 )  THEN
168          PRINT*, '+++ check_parameters:'
169          PRINT*, '    psolver="', TRIM( psolver ), '" only works for a ', &
170                       '1d domain-decomposition along x'
171          PRINT*, '    please do not set npey/=1 in the parameter file'
172       ENDIF
173       CALL local_stop
174    ENDIF
175    IF ( ( psolver == 'poisfft_hybrid'  .OR.  psolver == 'multigrid' )  .AND.  &
176         ( nxra > nxr  .OR.  nyna > nyn  .OR.  nza > nz ) )  THEN
177       IF ( myid == 0 )  THEN
178          PRINT*, '+++ check_parameters:'
179          PRINT*, '    psolver="', TRIM( psolver ), '" does not work for ', &
180                       'subdomains with unequal size'
181          PRINT*, '    please set grid_matching = ''strict'' in the parameter',&
182                       ' file'
183       ENDIF
184       CALL local_stop
185    ENDIF
186#else
187    IF ( psolver == 'poisfft_hybrid' )  THEN
188       IF ( myid == 0 )  THEN
189          PRINT*, '+++ check_parameters:'
190          PRINT*, '    psolver="', TRIM( psolver ), '" only works for a ', &
191                       'parallel environment'
192       ENDIF
193       CALL local_stop
194    ENDIF
195#endif
196
197    IF ( psolver == 'multigrid' )  THEN
198       IF ( cycle_mg == 'w' )  THEN
199          gamma_mg = 2
200       ELSEIF ( cycle_mg == 'v' )  THEN
201          gamma_mg = 1
202       ELSE
203          IF ( myid == 0 )  THEN
204             PRINT*, '+++ check_parameters:'
205             PRINT*, '    unknown multigrid cycle: cycle_mg=', cycle_mg
206          ENDIF
207          CALL local_stop
208       ENDIF
209    ENDIF
210
211    IF ( fft_method /= 'singleton-algorithm'  .AND.  &
212         fft_method /= 'temperton-algorithm'  .AND.  &
213         fft_method /= 'system-specific' )  THEN
214       IF ( myid == 0 )  THEN
215          PRINT*, '+++ check_parameters:'
216          PRINT*, '    unknown fft-algorithm: fft_method=', fft_method
217       ENDIF
218       CALL local_stop
219    ENDIF
220
221!
222!-- Advection schemes:
223    IF ( momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ups-scheme' ) &
224    THEN
225       IF ( myid == 0 )  PRINT*, '+++ check_parameters:  unknown advection ', &
226                                 'scheme: momentum_advec=', momentum_advec
227       CALL local_stop
228    ENDIF
229    IF ( ( momentum_advec == 'ups-scheme'  .OR.  scalar_advec == 'ups-scheme' )&
230                                      .AND.  timestep_scheme /= 'euler' )  THEN
231       IF ( myid == 0 )  PRINT*, '+++ check_parameters:  momentum_advec=', &
232                                 momentum_advec, ' is not allowed with ', &
233                                 'timestep_scheme=', timestep_scheme
234       CALL local_stop
235    ENDIF
236
237    IF ( scalar_advec /= 'pw-scheme'  .AND.  scalar_advec /= 'bc-scheme'  .AND.&
238         scalar_advec /= 'ups-scheme' )  THEN
239       IF ( myid == 0 )  PRINT*, '+++ check_parameters:  unknown advection ', &
240                                 'scheme: scalar_advec=', scalar_advec
241       CALL local_stop
242    ENDIF
243
244    IF ( use_sgs_for_particles  .AND.  .NOT. use_upstream_for_tke )  THEN
245       use_upstream_for_tke = .TRUE.
246       IF ( myid == 0 )  THEN
247          PRINT*, '+++ WARNING: check_parameters:  use_upstream_for_tke set ', &
248                       '.TRUE. because use_sgs_for_particles = .TRUE.'
249       ENDIF
250    ENDIF
251
252    IF ( use_upstream_for_tke  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
253       IF ( myid == 0 )  THEN
254          PRINT*, '+++ check_parameters:  use_upstream_for_tke = .TRUE. ', &
255                       'not allowed with timestep_scheme=', timestep_scheme
256       ENDIF
257       CALL local_stop
258    ENDIF
259
260!
261!-- Timestep schemes:
262    SELECT CASE ( TRIM( timestep_scheme ) )
263
264       CASE ( 'euler' )
265          intermediate_timestep_count_max = 1
266          asselin_filter_factor           = 0.0
267
268       CASE ( 'leapfrog', 'leapfrog+euler' )
269          intermediate_timestep_count_max = 1
270
271       CASE ( 'runge-kutta-2' )
272          intermediate_timestep_count_max = 2
273          asselin_filter_factor           = 0.0
274
275       CASE ( 'runge-kutta-3' )
276          intermediate_timestep_count_max = 3
277          asselin_filter_factor           = 0.0
278
279       CASE DEFAULT
280          IF ( myid == 0 )  PRINT*, '+++ check_parameters:  unknown timestep ',&
281                                    'scheme: timestep_scheme=', timestep_scheme
282          CALL local_stop
283
284    END SELECT
285
[63]286    IF ( scalar_advec == 'ups-scheme'  .AND.  timestep_scheme(1:5) == 'runge' )&
[1]287    THEN
288       IF ( myid == 0 )  THEN
289          PRINT*, '+++ check_parameters:  scalar advection scheme "', &
290                                          TRIM( scalar_advec ), '"'
291          PRINT*, '    does not work with timestep_scheme "', &
292                                          TRIM( timestep_scheme ), '"'
293       ENDIF
294       CALL local_stop
295    ENDIF
296
297    IF ( momentum_advec /= 'pw-scheme' .AND. timestep_scheme(1:5) == 'runge' ) &
298    THEN
299       IF ( myid == 0 )  THEN
300          PRINT*, '+++ check_parameters:  momentum advection scheme "', &
301                                          TRIM( momentum_advec ), '"'
302          PRINT*, '    does not work with timestep_scheme "', &
303                                          TRIM( timestep_scheme ), '"'
304       ENDIF
305       CALL local_stop
306    ENDIF
307
308
309    IF ( initializing_actions == ' ' )  THEN
310       IF ( myid == 0 )  THEN
311          PRINT*, '+++ check parameters:'
312          PRINT*, '    no value found for initializing_actions'
313       ENDIF
314       CALL local_stop
315    ENDIF
316
317    IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
318!
319!--    No model continuation run; several initialising actions are possible
320       action = initializing_actions
321       DO WHILE ( TRIM( action ) /= '' )
322          position = INDEX( action, ' ' )
323          SELECT CASE ( action(1:position-1) )
324
325             CASE ( 'set_constant_profiles', 'set_1d-model_profiles', &
[46]326                    'by_user', 'initialize_vortex',     'initialize_ptanom' )
[1]327                action = action(position+1:)
328
329             CASE DEFAULT
330                IF ( myid == 0 )  PRINT*, '+++ check_parameters: initializi', &
331                               'ng_action unkown or not allowed: action = "', &
332                               TRIM(action), '"'
333                CALL local_stop
334
335          END SELECT
336       ENDDO
337    ENDIF
338    IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0  .AND. &
339         INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
340       IF ( myid == 0 )  PRINT*, '+++ check_parameters: initializing_actions', &
341          '"set_constant_profiles" and "set_1d-model_profiles" are not', &
342          ' allowed simultaneously'
343       CALL local_stop
344    ENDIF
[46]345    IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0  .AND. &
346         INDEX( initializing_actions, 'by_user' ) /= 0 )  THEN
347       IF ( myid == 0 )  PRINT*, '+++ check_parameters: initializing_actions', &
348          '"set_constant_profiles" and "by_user" are not', &
349          ' allowed simultaneously'
350       CALL local_stop
351    ENDIF
352    IF ( INDEX( initializing_actions, 'by_user' ) /= 0  .AND. &
353         INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
354       IF ( myid == 0 )  PRINT*, '+++ check_parameters: initializing_actions', &
355          '"by_user" and "set_1d-model_profiles" are not', &
356          ' allowed simultaneously'
357       CALL local_stop
358    ENDIF
[1]359
[75]360    IF ( cloud_physics  .AND.  .NOT. humidity )  THEN
[72]361       IF ( myid == 0 )  PRINT*, '+++ check_parameters: cloud_physics =', &
362                                 cloud_physics, ' is not allowed with ',  &
[75]363                                 'humidity =', humidity
[1]364       CALL local_stop
365    ENDIF
366
[72]367    IF ( precipitation  .AND.  .NOT.  cloud_physics )  THEN
368       IF ( myid == 0 )  PRINT*, '+++ check_parameters: precipitation =', &
369                                 precipitation, ' is not allowed with ',  &
370                                 'cloud_physics =', cloud_physics
371       CALL local_stop
372    ENDIF
373
[75]374    IF ( humidity  .AND.  sloping_surface )  THEN
375       IF ( myid == 0 )  PRINT*, '+++ check_parameters: humidity = TRUE', &
[1]376                                 'and hang = TRUE are not',               &
377                                 ' allowed simultaneously' 
378       CALL local_stop       
379    ENDIF
380
[75]381    IF ( humidity  .AND.  scalar_advec == 'ups-scheme' )  THEN
[1]382       IF ( myid == 0 )  PRINT*, '+++ check_parameters: UPS-scheme', &
[75]383                                 'is not implemented for humidity'
[1]384       CALL local_stop       
385    ENDIF
386
[75]387    IF ( passive_scalar  .AND.  humidity )  THEN
388       IF ( myid == 0 )  PRINT*, '+++ check_parameters: humidity = TRUE and', &
[1]389                                 'passive_scalar = TRUE is not allowed ',     &
390                                 'simultaneously'
391       CALL local_stop
392    ENDIF
393
394    IF ( passive_scalar  .AND.  scalar_advec == 'ups-scheme' )  THEN
395       IF ( myid == 0 )  PRINT*, '+++ check_parameters: UPS-scheme', &
396                                 'is not implemented for passive_scalar'
397       CALL local_stop       
398    ENDIF
399
400    IF ( grid_matching /= 'strict'  .AND.  grid_matching /= 'match' )  THEN
401       IF ( myid == 0 )  PRINT*, '+++ check_parameters: illegal value "', &
402                                 TRIM( grid_matching ),                   &
403                                 '" found for parameter grid_matching'
404       CALL local_stop       
405    ENDIF
406
407!
408!-- In case of no model continuation run, check initialising parameters and
409!-- deduce further quantities
410    IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
411
412!
413!--    Initial profiles for 1D and 3D model, respectively
414       u_init  = ug_surface
415       v_init  = vg_surface
416       pt_init = pt_surface
[75]417       IF ( humidity )        q_init = q_surface
[1]418       IF ( passive_scalar )  q_init = s_surface
419
420!
421!--
422!--    If required, compute initial profile of the geostrophic wind
423!--    (component ug)
424       i = 1
425       gradient = 0.0
426       ug_vertical_gradient_level_ind(1) = 0
427       ug(0) = ug_surface
428       DO  k = 1, nzt+1
429          IF ( ug_vertical_gradient_level(i) < zu(k)  .AND. &
430               ug_vertical_gradient_level(i) >= 0.0 )  THEN
431             gradient = ug_vertical_gradient(i) / 100.0
432             ug_vertical_gradient_level_ind(i) = k - 1
433             i = i + 1
434             IF ( i > 10 )  THEN
435                IF ( myid == 0 )  THEN
436                   PRINT*, '+++ check_parameters: upper bound 10 of array', &
437                           ' "ug_vertical_gradient_level_ind" exceeded'
438                ENDIF
439                CALL local_stop
440             ENDIF
441          ENDIF
442          IF ( gradient /= 0.0 )  THEN
443             IF ( k /= 1 )  THEN
444                ug(k) = ug(k-1) + dzu(k) * gradient
445             ELSE
446                ug(k) = ug_surface + 0.5 * dzu(k) * gradient
447             ENDIF
448          ELSE
449             ug(k) = ug(k-1)
450          ENDIF
451       ENDDO
452
453       u_init = ug
454
455!
456!--    In case of no given gradients for ug, choose a vanishing gradient
457       IF ( ug_vertical_gradient_level(1) == -1.0 )  THEN
458          ug_vertical_gradient_level(1) = 0.0
459       ENDIF 
460
461!
462!--
463!--    If required, compute initial profile of the geostrophic wind
464!--    (component vg)
465       i = 1
466       gradient = 0.0
467       vg_vertical_gradient_level_ind(1) = 0
468       vg(0) = vg_surface
469       DO  k = 1, nzt+1
470          IF ( vg_vertical_gradient_level(i) < zu(k)  .AND. &
471               vg_vertical_gradient_level(i) >= 0.0 )  THEN
472             gradient = vg_vertical_gradient(i) / 100.0
473             vg_vertical_gradient_level_ind(i) = k - 1
474             i = i + 1
475             IF ( i > 10 )  THEN
476                IF ( myid == 0 )  THEN
477                   PRINT*, '+++ check_parameters: upper bound 10 of array', &
478                           ' "vg_vertical_gradient_level_ind" exceeded'
479                ENDIF
480                CALL local_stop
481             ENDIF
482          ENDIF
483          IF ( gradient /= 0.0 )  THEN
484             IF ( k /= 1 )  THEN
485                vg(k) = vg(k-1) + dzu(k) * gradient
486             ELSE
487                vg(k) = vg_surface + 0.5 * dzu(k) * gradient
488             ENDIF
489          ELSE
490             vg(k) = vg(k-1)
491          ENDIF
492       ENDDO
493
494       v_init = vg
495 
496!
497!--    In case of no given gradients for vg, choose a vanishing gradient
498       IF ( vg_vertical_gradient_level(1) == -1.0 )  THEN
499          vg_vertical_gradient_level(1) = 0.0
500       ENDIF
501
502!
503!--    If required, compute initial temperature profile using the given
504!--    temperature gradient
505       i = 1
506       gradient = 0.0
507       pt_vertical_gradient_level_ind(1) = 0
508       DO  k = 1, nzt+1
509          IF ( pt_vertical_gradient_level(i) < zu(k)  .AND. &
510               pt_vertical_gradient_level(i) >= 0.0 )  THEN
511             gradient = pt_vertical_gradient(i) / 100.0
512             pt_vertical_gradient_level_ind(i) = k - 1
513             i = i + 1
514             IF ( i > 10 )  THEN
515                IF ( myid == 0 )  THEN
516                   PRINT*, '+++ check_parameters: upper bound 10 of array', &
517                           ' "pt_vertical_gradient_level_ind" exceeded'
518                ENDIF
519                CALL local_stop
520             ENDIF
521          ENDIF
522          IF ( gradient /= 0.0 )  THEN
523             IF ( k /= 1 )  THEN
524                pt_init(k) = pt_init(k-1) + dzu(k) * gradient
525             ELSE
526                pt_init(k) = pt_init(k-1) + 0.5 * dzu(k) * gradient
527             ENDIF
528          ELSE
529             pt_init(k) = pt_init(k-1)
530          ENDIF
531       ENDDO
532
533!
534!--    In case of no given temperature gradients, choose gradient of neutral
535!--    stratification
536       IF ( pt_vertical_gradient_level(1) == -1.0 )  THEN
537          pt_vertical_gradient_level(1) = 0.0
538       ENDIF
539
540!
541!--    Store temperature gradient at the top boundary for possile Neumann
542!--    boundary condition
[19]543       bc_pt_t_val = ( pt_init(nzt+1) - pt_init(nzt) ) / dzu(nzt+1)
[1]544
545!
546!--    If required, compute initial humidity or scalar profile using the given
547!--    humidity/scalar gradient. In case of scalar transport, initially store
548!--    values of the scalar parameters on humidity parameters
549       IF ( passive_scalar )  THEN
550          bc_q_b                    = bc_s_b
551          bc_q_t                    = bc_s_t
552          q_surface                 = s_surface
553          q_surface_initial_change  = s_surface_initial_change
554          q_vertical_gradient       = s_vertical_gradient
555          q_vertical_gradient_level = s_vertical_gradient_level
556          surface_waterflux         = surface_scalarflux
557       ENDIF
558
[75]559       IF ( humidity  .OR.  passive_scalar )  THEN
[1]560
561          i = 1
562          gradient = 0.0
563          q_vertical_gradient_level_ind(1) = 0
564          DO  k = 1, nzt+1
565             IF ( q_vertical_gradient_level(i) < zu(k)  .AND. &
566                  q_vertical_gradient_level(i) >= 0.0 )  THEN
567                gradient = q_vertical_gradient(i) / 100.0
568                q_vertical_gradient_level_ind(i) = k - 1
569                i = i + 1
570                IF ( i > 10 )  THEN
571                   IF ( myid == 0 )  THEN
572                      PRINT*, '+++ check_parameters: upper bound 10 of arr', &
573                              'ay "q_vertical_gradient_level_ind" exceeded'
574                   ENDIF
575                   CALL local_stop
576                ENDIF
577             ENDIF
578             IF ( gradient /= 0.0 )  THEN
579                IF ( k /= 1 )  THEN
580                   q_init(k) = q_init(k-1) + dzu(k) * gradient
581                ELSE
582                   q_init(k) = q_init(k-1) + 0.5 * dzu(k) * gradient
583                ENDIF
584             ELSE
585                q_init(k) = q_init(k-1)
586             ENDIF
[72]587!
588!--          Avoid negative humidities
589             IF ( q_init(k) < 0.0 )  THEN
590                q_init(k) = 0.0
591             ENDIF
[1]592          ENDDO
593
594!
595!--       In case of no given humidity gradients, choose zero gradient
596!--       conditions
597          IF ( q_vertical_gradient_level(1) == -1.0 )  THEN
598             q_vertical_gradient_level(1) = 0.0
599          ENDIF
600
601!
602!--       Store humidity gradient at the top boundary for possile Neumann
603!--       boundary condition
[19]604          bc_q_t_val = ( q_init(nzt+1) - q_init(nzt) ) / dzu(nzt+1)
[1]605
606       ENDIF
607
608    ENDIF
609
610!
611!-- Compute Coriolis parameter
612    f  = 2.0 * omega * SIN( phi / 180.0 * pi )
613    fs = 2.0 * omega * COS( phi / 180.0 * pi )
614
615!
[57]616!-- Reference temperature to be used in buoyancy terms
617    IF ( pt_reference /= 9999999.9 )  use_pt_reference = .TRUE.
618
619!
[1]620!-- In case of a given slope, compute the relevant quantities
621    IF ( alpha_surface /= 0.0 )  THEN
622       IF ( ABS( alpha_surface ) > 90.0 )  THEN
623          IF ( myid == 0 )  PRINT*, '+++ check_parameters: ABS( alpha_surface',&
624                                    '=', alpha_surface, ' ) must be < 90.0'
625          CALL local_stop
626       ENDIF
627       sloping_surface = .TRUE.
628       cos_alpha_surface = COS( alpha_surface / 180.0 * pi )
629       sin_alpha_surface = SIN( alpha_surface / 180.0 * pi )
630    ENDIF
631
632!
633!-- Check time step and cfl_factor
634    IF ( dt /= -1.0 )  THEN
635       IF ( dt <= 0.0  .AND.  dt /= -1.0 )  THEN
636          IF ( myid == 0 )  PRINT*, '+++ check_parameters:  dt=', dt, ' <= 0.0'
637          CALL local_stop
638       ENDIF
639       dt_3d = dt
640       dt_fixed = .TRUE.
641    ENDIF
642
643    IF ( cfl_factor <= 0.0  .OR.  cfl_factor > 1.0 )  THEN
644       IF ( cfl_factor == -1.0 )  THEN
645          IF ( momentum_advec == 'ups-scheme'  .OR.  &
646               scalar_advec == 'ups-scheme' )  THEN
647             cfl_factor = 0.8
648          ELSE
649             IF ( timestep_scheme == 'runge-kutta-2' )  THEN
650                cfl_factor = 0.8
651             ELSEIF ( timestep_scheme == 'runge-kutta-3' )  THEN
652                cfl_factor = 0.9
653             ELSE
654                cfl_factor = 0.1
655             ENDIF
656          ENDIF
657       ELSE
658          IF ( myid == 0 )  THEN
659             PRINT*, '+++ check_parameters: cfl_factor=', cfl_factor, &
660                         ' out of range'
661             PRINT*, '+++                   0.0 < cfl_factor <= 1.0 is required'
662          ENDIF
663          CALL local_stop
664       ENDIF
665    ENDIF
666
667!
668!-- Store simulated time at begin
669    simulated_time_at_begin = simulated_time
670
671!
672!-- Set wind speed in the Galilei-transformed system
673    IF ( galilei_transformation )  THEN
674       IF ( use_ug_for_galilei_tr .AND.                &
675            ug_vertical_gradient_level(1) == 0.0 .AND. & 
676            vg_vertical_gradient_level(1) == 0.0 )  THEN
677          u_gtrans = ug_surface
678          v_gtrans = vg_surface
679       ELSEIF ( use_ug_for_galilei_tr .AND.                &
680                ug_vertical_gradient_level(1) /= 0.0 )  THEN
681          IF ( myid == 0 )  THEN
682             PRINT*, '+++ check_parameters:'
683             PRINT*, '    baroclinicity (ug) not allowed'
684             PRINT*, '    simultaneously with galilei transformation'
685          ENDIF
686          CALL local_stop
687       ELSEIF ( use_ug_for_galilei_tr .AND.                &
688                vg_vertical_gradient_level(1) /= 0.0 )  THEN
689          IF ( myid == 0 )  THEN
690             PRINT*, '+++ check_parameters:'
691             PRINT*, '    baroclinicity (vg) not allowed'
692             PRINT*, '    simultaneously with galilei transformation'
693          ENDIF
694          CALL local_stop
695       ELSE
696          IF ( myid == 0 )  THEN
697             PRINT*, '+++ WARNING: check_parameters:'
698             PRINT*, '    variable translation speed used for galilei-tran' // &
699                          'sformation, which'
700             PRINT*, '    may cause instabilities in stably stratified regions'
701          ENDIF
702       ENDIF
703    ENDIF
704
705!
706!-- In case of using a prandtl-layer, calculated (or prescribed) surface
707!-- fluxes have to be used in the diffusion-terms
708    IF ( prandtl_layer )  use_surface_fluxes = .TRUE.
709
710!
711!-- Check boundary conditions and set internal variables:
712!-- Lateral boundary conditions
[73]713    IF ( bc_lr /= 'cyclic'  .AND.  bc_lr /= 'dirichlet/radiation'  .AND. &
714         bc_lr /= 'radiation/dirichlet' )  THEN
[1]715       IF ( myid == 0 )  THEN
716          PRINT*, '+++ check_parameters:'
717          PRINT*, '    unknown boundary condition: bc_lr = ', bc_lr
718       ENDIF
719       CALL local_stop
720    ENDIF
[73]721    IF ( bc_ns /= 'cyclic'  .AND.  bc_ns /= 'dirichlet/radiation'  .AND. &
722         bc_ns /= 'radiation/dirichlet' )  THEN
[1]723       IF ( myid == 0 )  THEN
724          PRINT*, '+++ check_parameters:'
725          PRINT*, '    unknown boundary condition: bc_ns = ', bc_ns
726       ENDIF
727       CALL local_stop
728    ENDIF
729
730!
731!-- Non-cyclic lateral boundaries require the multigrid method and Piascek-
732!-- Willimas advection scheme. Several schemes and tools do not work with
733!-- non-cyclic boundary conditions.
734    IF ( bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic' )  THEN
735       IF ( psolver /= 'multigrid' )  THEN
736          IF ( myid == 0 )  THEN
737             PRINT*, '+++ check_parameters:'
738             PRINT*, '    non-cyclic lateral boundaries do not allow', &
739                          ' psolver = ', psolver 
740          ENDIF
741          CALL local_stop
742       ENDIF
743       IF ( momentum_advec /= 'pw-scheme' )  THEN
744          IF ( myid == 0 )  THEN
745             PRINT*, '+++ check_parameters:'
746             PRINT*, '    non-cyclic lateral boundaries do not allow', &
747                          ' momentum_advec = ', momentum_advec 
748          ENDIF
749          CALL local_stop
750       ENDIF
751       IF ( scalar_advec /= 'pw-scheme' )  THEN
752          IF ( myid == 0 )  THEN
753             PRINT*, '+++ check_parameters:'
754             PRINT*, '    non-cyclic lateral boundaries do not allow', &
755                          ' scalar_advec = ', scalar_advec 
756          ENDIF
757          CALL local_stop
758       ENDIF
759       IF ( galilei_transformation )  THEN
760          IF ( myid == 0 )  THEN
761             PRINT*, '+++ check_parameters:'
762             PRINT*, '    non-cyclic lateral boundaries do not allow', &
763                          ' galilei_transformation = .T.' 
764          ENDIF
765          CALL local_stop
766       ENDIF
[75]767!       IF ( conserve_volume_flow )  THEN
768!          IF ( myid == 0 )  THEN
769!             PRINT*, '+++ check_parameters:'
770!             PRINT*, '    non-cyclic lateral boundaries do not allow', &
771!                          ' conserve_volume_flow = .T.'
772!          ENDIF
773!          CALL local_stop
774!       ENDIF
[1]775    ENDIF
776
777!
778!-- Bottom boundary condition for the turbulent Kinetic energy
779    IF ( bc_e_b == 'neumann' )  THEN
780       ibc_e_b = 1
781       IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
782          IF ( myid == 0 )  THEN
783             PRINT*, '+++ WARNING: check_parameters:'
784             PRINT*, '    adjust_mixing_length = TRUE and bc_e_b = ', bc_e_b
785          ENDIF
786       ENDIF
787    ELSEIF ( bc_e_b == '(u*)**2+neumann' )  THEN
788       ibc_e_b = 2
789       IF ( .NOT. adjust_mixing_length  .AND.  prandtl_layer )  THEN
790          IF ( myid == 0 )  THEN
791             PRINT*, '+++ WARNING: check_parameters:'
792             PRINT*, '    adjust_mixing_length = FALSE and bc_e_b = ', bc_e_b
793          ENDIF
794       ENDIF
795       IF ( .NOT. prandtl_layer )  THEN
796          bc_e_b = 'neumann'
797          ibc_e_b = 1
798          IF ( myid == 0 )  THEN
799             PRINT*, '+++ WARNING: check_parameters:'
800             PRINT*, '    boundary condition bc_e_b changed to "', bc_e_b, '"'
801          ENDIF
802       ENDIF
803    ELSE
804       IF ( myid == 0 )  THEN
805          PRINT*, '+++ check_parameters:'
806          PRINT*, '    unknown boundary condition: bc_e_b = ', bc_e_b
807       ENDIF
808       CALL local_stop
809    ENDIF
810
811!
812!-- Boundary conditions for perturbation pressure
813    IF ( bc_p_b == 'dirichlet' )  THEN
814       ibc_p_b = 0
815    ELSEIF ( bc_p_b == 'neumann' )  THEN
816       ibc_p_b = 1
817    ELSEIF ( bc_p_b == 'neumann+inhomo' )  THEN
818       ibc_p_b = 2
819    ELSE
820       IF ( myid == 0 )  THEN
821          PRINT*, '+++ check_parameters:'
822          PRINT*, '    unknown boundary condition: bc_p_b = ', bc_p_b
823       ENDIF
824       CALL local_stop
825    ENDIF
826    IF ( ibc_p_b == 2  .AND.  .NOT. prandtl_layer )  THEN
827       IF ( myid == 0 )  THEN
828          PRINT*, '+++ check_parameters:'
829          PRINT*, '    boundary condition: bc_p_b = ', TRIM( bc_p_b ), &
830                       ' not allowed with'
831          PRINT*, '    prandtl_layer = .FALSE.' 
832       ENDIF
833       CALL local_stop
834    ENDIF
835    IF ( bc_p_t == 'dirichlet' )  THEN
836       ibc_p_t = 0
837    ELSEIF ( bc_p_t == 'neumann' )  THEN
838       ibc_p_t = 1
839    ELSE
840       IF ( myid == 0 )  THEN
841          PRINT*, '+++ check_parameters:'
842          PRINT*, '    unknown boundary condition: bc_p_t = ', bc_p_t
843       ENDIF
844       CALL local_stop
845    ENDIF
846
847!
848!-- Boundary conditions for potential temperature
849    IF ( bc_pt_b == 'dirichlet' )  THEN
850       ibc_pt_b = 0
851    ELSEIF ( bc_pt_b == 'neumann' )  THEN
852       ibc_pt_b = 1
853    ELSE
854       IF ( myid == 0 )  THEN
855          PRINT*, '+++ check_parameters:'
856          PRINT*, '    unknown boundary condition: bc_pt_b = ', bc_pt_b
857       ENDIF
858       CALL local_stop
859    ENDIF
860    IF ( bc_pt_t == 'dirichlet' )  THEN
861       ibc_pt_t = 0
862    ELSEIF ( bc_pt_t == 'neumann' )  THEN
863       ibc_pt_t = 1
[19]864    ELSEIF ( bc_pt_t == 'initial_gradient' )  THEN
865       ibc_pt_t = 2
[1]866    ELSE
867       IF ( myid == 0 )  THEN
868          PRINT*, '+++ check_parameters:'
869          PRINT*, '    unknown boundary condition: bc_pt_t = ', bc_pt_t
870       ENDIF
871       CALL local_stop
872    ENDIF
873
[20]874    IF ( surface_heatflux == 9999999.9 )  constant_heatflux     = .FALSE.
875    IF ( top_heatflux     == 9999999.9 )  constant_top_heatflux = .FALSE.
[1]876
877!
878!-- A given surface temperature implies Dirichlet boundary condition for
879!-- temperature. In this case specification of a constant heat flux is
880!-- forbidden.
881    IF ( ibc_pt_b == 0  .AND.   constant_heatflux  .AND. &
882         surface_heatflux /= 0.0 )  THEN
883       IF ( myid == 0 )  THEN
884          PRINT*, '+++ check_parameters:'
885          PRINT*, '    boundary_condition: bc_pt_b = ', bc_pt_b
886          PRINT*, '    is not allowed with constant_heatflux = .TRUE.'
887       ENDIF
888       CALL local_stop
889    ENDIF
890    IF ( constant_heatflux  .AND.  pt_surface_initial_change /= 0.0 )  THEN
891       IF ( myid == 0 )  THEN
892          PRINT*, '+++ check_parameters: constant_heatflux = .TRUE. is not'
893          PRINT*, '    allowed with pt_surface_initial_change (/=0) = ', &
894                  pt_surface_initial_change
895       ENDIF
896       CALL local_stop
897    ENDIF
898
899!
[19]900!-- A given temperature at the top implies Dirichlet boundary condition for
901!-- temperature. In this case specification of a constant heat flux is
902!-- forbidden.
903    IF ( ibc_pt_t == 0  .AND.   constant_top_heatflux  .AND. &
904         top_heatflux /= 0.0 )  THEN
905       IF ( myid == 0 )  THEN
906          PRINT*, '+++ check_parameters:'
907          PRINT*, '    boundary_condition: bc_pt_t = ', bc_pt_t
908          PRINT*, '    is not allowed with constant_top_heatflux = .TRUE.'
909       ENDIF
910       CALL local_stop
911    ENDIF
912
913!
[75]914!-- In case of humidity or passive scalar, set boundary conditions for total
[1]915!-- water content / scalar
[75]916    IF ( humidity  .OR.  passive_scalar ) THEN
917       IF ( humidity )  THEN
[1]918          sq = 'q'
919       ELSE
920          sq = 's'
921       ENDIF
922       IF ( bc_q_b == 'dirichlet' )  THEN
923          ibc_q_b = 0
924       ELSEIF ( bc_q_b == 'neumann' )  THEN
925          ibc_q_b = 1
926       ELSE
927          IF ( myid == 0 )  THEN
928             PRINT*, '+++ check_parameters:'
929             PRINT*, '    unknown boundary condition: bc_', sq, '_b = ', bc_q_b
930          ENDIF
931          CALL local_stop
932       ENDIF
933       IF ( bc_q_t == 'dirichlet' )  THEN
934          ibc_q_t = 0
935       ELSEIF ( bc_q_t == 'neumann' )  THEN
936          ibc_q_t = 1
937       ELSE
938          IF ( myid == 0 )  THEN
939             PRINT*, '+++ check_parameters:'
940             PRINT*, '    unknown boundary condition: bc_', sq, '_t = ', bc_q_t
941          ENDIF
942          CALL local_stop
943       ENDIF
944
945       IF ( surface_waterflux == 0.0 )  constant_waterflux = .FALSE.
946
947!
948!--    A given surface humidity implies Dirichlet boundary condition for
[75]949!--    humidity. In this case specification of a constant water flux is
[1]950!--    forbidden.
951       IF ( ibc_q_b == 0  .AND.  constant_waterflux )  THEN
952          IF ( myid == 0 )  THEN
953             PRINT*, '+++ check_parameters:'
954             PRINT*, '    boundary_condition: bc_', sq, '_b = ', bc_q_b
955             PRINT*, '    is not allowed with prescribed surface flux'
956          ENDIF
957          CALL local_stop
958       ENDIF
959       IF ( constant_waterflux  .AND.  q_surface_initial_change /= 0.0 )  THEN
960          IF ( myid == 0 )  THEN
961             PRINT*, '+++ check_parameters: a prescribed surface flux is not'
962             PRINT*, '    allowed with ', sq, '_surface_initial_change (/=0)', &
963                     ' = ', q_surface_initial_change
964          ENDIF
965          CALL local_stop
966       ENDIF
967       
968    ENDIF
969
970!
971!-- Boundary conditions for horizontal components of wind speed
972    IF ( bc_uv_b == 'dirichlet' )  THEN
973       ibc_uv_b = 0
974    ELSEIF ( bc_uv_b == 'neumann' )  THEN
975       ibc_uv_b = 1
976       IF ( prandtl_layer )  THEN
977          IF ( myid == 0 )  THEN
978             PRINT*, '+++ check_parameters:'
979             PRINT*, '    boundary condition: bc_uv_b = ', TRIM( bc_uv_b ), &
980                          ' is not allowed with'
981             PRINT*, '    prandtl_layer = .TRUE.' 
982          ENDIF
983          CALL local_stop
984       ENDIF
985    ELSE
986       IF ( myid == 0 )  THEN
987          PRINT*, '+++ check_parameters:'
988          PRINT*, '    unknown boundary condition: bc_uv_b = ', bc_uv_b
989       ENDIF
990       CALL local_stop
991    ENDIF
992    IF ( bc_uv_t == 'dirichlet' )  THEN
993       ibc_uv_t = 0
994    ELSEIF ( bc_uv_t == 'neumann' )  THEN
995       ibc_uv_t = 1
996    ELSE
997       IF ( myid == 0 )  THEN
998          PRINT*, '+++ check_parameters:'
999          PRINT*, '    unknown boundary condition: bc_uv_t = ', bc_uv_t
1000       ENDIF
1001       CALL local_stop
1002    ENDIF
1003
1004!
1005!-- Compute and check, respectively, the Rayleigh Damping parameter
1006    IF ( rayleigh_damping_factor == -1.0 )  THEN
1007       IF ( momentum_advec == 'ups-scheme' )  THEN
1008          rayleigh_damping_factor = 0.01
1009       ELSE
1010          rayleigh_damping_factor = 0.0
1011       ENDIF
1012    ELSE
1013       IF ( rayleigh_damping_factor < 0.0 .OR. rayleigh_damping_factor > 1.0 ) &
1014       THEN
1015          IF ( myid == 0 )  THEN
1016             PRINT*, '+++ check_parameters:'
1017             PRINT*, '    rayleigh_damping_factor = ', rayleigh_damping_factor,&
1018                          ' out of range [0.0,1.0]'
1019          ENDIF
1020          CALL local_stop
1021       ENDIF
1022    ENDIF
1023
1024    IF ( rayleigh_damping_height == -1.0 )  THEN
1025       rayleigh_damping_height = 0.66666666666 * zu(nzt)
1026    ELSE
1027       IF ( rayleigh_damping_height < 0.0  .OR. &
1028            rayleigh_damping_height > zu(nzt) )  THEN
1029          IF ( myid == 0 )  THEN
1030             PRINT*, '+++ check_parameters:'
1031             PRINT*, '    rayleigh_damping_height = ', rayleigh_damping_height,&
1032                          ' out of range [0.0,', zu(nzt), ']'
1033          ENDIF
1034          CALL local_stop
1035       ENDIF
1036    ENDIF
1037
1038!
1039!-- Check limiters for Upstream-Spline scheme
1040    IF ( overshoot_limit_u < 0.0  .OR.  overshoot_limit_v < 0.0  .OR.  &
1041         overshoot_limit_w < 0.0  .OR.  overshoot_limit_pt < 0.0  .OR. &
1042         overshoot_limit_e < 0.0 )  THEN
1043       IF ( myid == 0 )  THEN
1044          PRINT*, '+++ check_parameters:'
1045          PRINT*, '    overshoot_limit_... < 0.0 is not allowed'
1046       ENDIF
1047       CALL local_stop
1048    ENDIF
1049    IF ( ups_limit_u < 0.0 .OR. ups_limit_v < 0.0 .OR. ups_limit_w < 0.0 .OR. &
1050         ups_limit_pt < 0.0 .OR. ups_limit_e < 0.0 )  THEN
1051       IF ( myid == 0 )  THEN
1052          PRINT*, '+++ check_parameters:'
1053          PRINT*, '    ups_limit_... < 0.0 is not allowed'
1054       ENDIF
1055       CALL local_stop
1056    ENDIF
1057
1058!
1059!-- Check number of chosen statistic regions. More than 10 regions are not
1060!-- allowed, because so far no more than 10 corresponding output files can
1061!-- be opened (cf. check_open)
1062    IF ( statistic_regions > 9  .OR.  statistic_regions < 0 )  THEN
1063       IF ( myid == 0 )  THEN
1064          PRINT*, '+++ check_parameters: Number of statistic_regions = ', &
1065                       statistic_regions+1
1066          PRINT*, '    Only 10 regions are allowed'
1067       ENDIF
1068       CALL local_stop
1069    ENDIF
1070    IF ( normalizing_region > statistic_regions  .OR. &
1071         normalizing_region < 0)  THEN
1072       IF ( myid == 0 )  THEN
1073          PRINT*, '+++ check_parameters: normalizing_region = ', &
1074                       normalizing_region, ' is unknown'
1075          PRINT*, '    Must be <= ', statistic_regions
1076       ENDIF
1077       CALL local_stop
1078    ENDIF
1079
1080!
1081!-- Set the default intervals for data output, if necessary
1082!-- NOTE: dt_dosp has already been set in package_parin
1083    IF ( dt_data_output /= 9999999.9 )  THEN
1084       IF ( dt_dopr           == 9999999.9 )  dt_dopr           = dt_data_output
1085       IF ( dt_dopts          == 9999999.9 )  dt_dopts          = dt_data_output
1086       IF ( dt_do2d_xy        == 9999999.9 )  dt_do2d_xy        = dt_data_output
1087       IF ( dt_do2d_xz        == 9999999.9 )  dt_do2d_xz        = dt_data_output
1088       IF ( dt_do2d_yz        == 9999999.9 )  dt_do2d_yz        = dt_data_output
1089       IF ( dt_do3d           == 9999999.9 )  dt_do3d           = dt_data_output
1090       IF ( dt_data_output_av == 9999999.9 )  dt_data_output_av = dt_data_output
1091    ENDIF
1092
1093!
1094!-- Set the default skip time intervals for data output, if necessary
1095    IF ( skip_time_dopr    == 9999999.9 ) &
1096                                       skip_time_dopr    = skip_time_data_output
1097    IF ( skip_time_dosp    == 9999999.9 ) &
1098                                       skip_time_dosp    = skip_time_data_output
1099    IF ( skip_time_do2d_xy == 9999999.9 ) &
1100                                       skip_time_do2d_xy = skip_time_data_output
1101    IF ( skip_time_do2d_xz == 9999999.9 ) &
1102                                       skip_time_do2d_xz = skip_time_data_output
1103    IF ( skip_time_do2d_yz == 9999999.9 ) &
1104                                       skip_time_do2d_yz = skip_time_data_output
1105    IF ( skip_time_do3d    == 9999999.9 ) &
1106                                       skip_time_do3d    = skip_time_data_output
1107    IF ( skip_time_data_output_av == 9999999.9 ) &
1108                                skip_time_data_output_av = skip_time_data_output
1109
1110!
1111!-- Check the average intervals (first for 3d-data, then for profiles and
1112!-- spectra)
1113    IF ( averaging_interval > dt_data_output_av )  THEN
1114       IF ( myid == 0 )  THEN
1115          PRINT*, '+++ check_parameters: average_interval=',              &
1116                       averaging_interval, ' must be <= dt_data_output=', &
1117                       dt_data_output
1118       ENDIF
1119       CALL local_stop
1120    ENDIF
1121
1122    IF ( averaging_interval_pr == 9999999.9 )  THEN
1123       averaging_interval_pr = averaging_interval
1124    ENDIF
1125
1126    IF ( averaging_interval_pr > dt_dopr )  THEN
1127       IF ( myid == 0 )  THEN
1128          PRINT*, '+++ check_parameters: averaging_interval_pr=', &
1129                       averaging_interval_pr, ' must be <= dt_dopr=', dt_dopr
1130       ENDIF
1131       CALL local_stop
1132    ENDIF
1133
1134    IF ( averaging_interval_sp == 9999999.9 )  THEN
1135       averaging_interval_sp = averaging_interval
1136    ENDIF
1137
1138    IF ( averaging_interval_sp > dt_dosp )  THEN
1139       IF ( myid == 0 )  THEN
1140          PRINT*, '+++ check_parameters: averaging_interval_sp=', &
1141                       averaging_interval_sp, ' must be <= dt_dosp=', dt_dosp
1142       ENDIF
1143       CALL local_stop
1144    ENDIF
1145
1146!
1147!-- Set the default interval for profiles entering the temporal average
1148    IF ( dt_averaging_input_pr == 9999999.9 )  THEN
1149       dt_averaging_input_pr = dt_averaging_input
1150    ENDIF
1151
1152!
1153!-- Set the default interval for the output of timeseries to a reasonable
1154!-- value (tries to minimize the number of calls of flow_statistics)
1155    IF ( dt_dots == 9999999.9 )  THEN
1156       IF ( averaging_interval_pr == 0.0 )  THEN
1157          dt_dots = MIN( dt_run_control, dt_dopr )
1158       ELSE
1159          dt_dots = MIN( dt_run_control, dt_averaging_input_pr )
1160       ENDIF
1161    ENDIF
1162
1163!
1164!-- Check the sample rate for averaging (first for 3d-data, then for profiles)
1165    IF ( dt_averaging_input > averaging_interval )  THEN
1166       IF ( myid == 0 )  THEN
1167          PRINT*, '+++ check_parameters: dt_averaging_input=',                &
1168                       dt_averaging_input, ' must be <= averaging_interval=', &
1169                       averaging_interval
1170       ENDIF
1171       CALL local_stop
1172    ENDIF
1173
1174    IF ( dt_averaging_input_pr > averaging_interval_pr )  THEN
1175       IF ( myid == 0 )  THEN
1176          PRINT*, '+++ check_parameters: dt_averaging_input_pr=', &
1177                       dt_averaging_input_pr,                     &
1178                       ' must be <= averaging_interval_pr=',      &
1179                       averaging_interval_pr
1180       ENDIF
1181       CALL local_stop
1182    ENDIF
1183
1184!
[72]1185!-- Set the default value for the integration interval of precipitation amount
1186    IF ( precipitation )  THEN
1187       IF ( precipitation_amount_interval == 9999999.9 )  THEN
1188          precipitation_amount_interval = dt_do2d_xy
1189       ELSE
1190          IF ( precipitation_amount_interval > dt_do2d_xy )  THEN
1191             IF ( myid == 0 )  PRINT*, '+++ check_parameters: ',              &
1192                                       'precipitation_amount_interval =',     &
1193                                        precipitation_amount_interval,        &
1194                                       ' must not be larger than dt_do2d_xy', &
1195                                       ' = ', dt_do2d_xy   
1196       CALL local_stop
1197          ENDIF
1198       ENDIF
1199    ENDIF
1200
1201!
[1]1202!-- Determine the number of output profiles and check whether they are
1203!-- permissible
1204    DO  WHILE ( data_output_pr(dopr_n+1) /= '          ' )
1205
1206       dopr_n = dopr_n + 1
1207       i = dopr_n
1208
1209!
1210!--    Determine internal profile number (for hom, homs)
1211!--    and store height levels
1212       SELECT CASE ( TRIM( data_output_pr(i) ) )
1213
1214          CASE ( 'u', '#u' )
1215             dopr_index(i) = 1
[87]1216             dopr_unit(i)  = 'm/s'
[1]1217             hom(:,2,1,:)  = SPREAD( zu, 2, statistic_regions+1 )
1218             IF ( data_output_pr(i)(1:1) == '#' )  THEN
1219                dopr_initial_index(i) = 5
1220                hom(:,2,5,:)          = SPREAD( zu, 2, statistic_regions+1 )
1221                data_output_pr(i)     = data_output_pr(i)(2:)
1222             ENDIF
1223
1224          CASE ( 'v', '#v' )
1225             dopr_index(i) = 2
[87]1226             dopr_unit(i)  = 'm/s'
1227             hom(:,2,2,:)  = SPREAD( zu, 2, statistic_regions+1 )
[1]1228             IF ( data_output_pr(i)(1:1) == '#' )  THEN
1229                dopr_initial_index(i) = 6
1230                hom(:,2,6,:)          = SPREAD( zu, 2, statistic_regions+1 )
1231                data_output_pr(i)     = data_output_pr(i)(2:)
1232             ENDIF
1233
1234          CASE ( 'w' )
1235             dopr_index(i) = 3
[87]1236             dopr_unit(i)  = 'm/s'
1237             hom(:,2,3,:)  = SPREAD( zw, 2, statistic_regions+1 )
[1]1238
1239          CASE ( 'pt', '#pt' )
1240             IF ( .NOT. cloud_physics ) THEN
1241                dopr_index(i) = 4
[87]1242                dopr_unit(i)  = 'K'
[1]1243                hom(:,2,4,:)  = SPREAD( zu, 2, statistic_regions+1 )
1244                IF ( data_output_pr(i)(1:1) == '#' )  THEN
1245                   dopr_initial_index(i) = 7
1246                   hom(:,2,7,:)          = SPREAD( zu, 2, statistic_regions+1 )
[87]1247                   hom(nzb,2,7,:)        = 0.0    ! because zu(nzb) is negative
[1]1248                   data_output_pr(i)     = data_output_pr(i)(2:)
1249                ENDIF
1250             ELSE
1251                dopr_index(i) = 43
[87]1252                dopr_unit(i)  = 'K'
[1]1253                hom(:,2,43,:)  = SPREAD( zu, 2, statistic_regions+1 )
1254                IF ( data_output_pr(i)(1:1) == '#' )  THEN
1255                   dopr_initial_index(i) = 28
1256                   hom(:,2,28,:)         = SPREAD( zu, 2, statistic_regions+1 )
[87]1257                   hom(nzb,2,28,:)       = 0.0    ! because zu(nzb) is negative
[1]1258                   data_output_pr(i)     = data_output_pr(i)(2:)
1259                ENDIF
1260             ENDIF
1261
1262          CASE ( 'e' )
1263             dopr_index(i)  = 8
[87]1264             dopr_unit(i)   = 'm2/s2'
[1]1265             hom(:,2,8,:)   = SPREAD( zu, 2, statistic_regions+1 )
1266             hom(nzb,2,8,:) = 0.0
1267
1268          CASE ( 'km', '#km' )
1269             dopr_index(i)  = 9
[87]1270             dopr_unit(i)   = 'm2/s'
[1]1271             hom(:,2,9,:)   = SPREAD( zu, 2, statistic_regions+1 )
1272             hom(nzb,2,9,:) = 0.0
1273             IF ( data_output_pr(i)(1:1) == '#' )  THEN
1274                dopr_initial_index(i) = 23
1275                hom(:,2,23,:)         = hom(:,2,9,:)
1276                data_output_pr(i)     = data_output_pr(i)(2:)
1277             ENDIF
1278
1279          CASE ( 'kh', '#kh' )
1280             dopr_index(i)   = 10
[87]1281             dopr_unit(i)    = 'm2/s'
[1]1282             hom(:,2,10,:)   = SPREAD( zu, 2, statistic_regions+1 )
1283             hom(nzb,2,10,:) = 0.0
1284             IF ( data_output_pr(i)(1:1) == '#' )  THEN
1285                dopr_initial_index(i) = 24
1286                hom(:,2,24,:)         = hom(:,2,10,:)
1287                data_output_pr(i)     = data_output_pr(i)(2:)
1288             ENDIF
1289
1290          CASE ( 'l', '#l' )
1291             dopr_index(i)   = 11
[87]1292             dopr_unit(i)    = 'm'
[1]1293             hom(:,2,11,:)   = SPREAD( zu, 2, statistic_regions+1 )
1294             hom(nzb,2,11,:) = 0.0
1295             IF ( data_output_pr(i)(1:1) == '#' )  THEN
1296                dopr_initial_index(i) = 25
1297                hom(:,2,25,:)         = hom(:,2,11,:)
1298                data_output_pr(i)     = data_output_pr(i)(2:)
1299             ENDIF
1300
1301          CASE ( 'w"u"' )
1302             dopr_index(i) = 12
[87]1303             dopr_unit(i)  = 'm2/s2'
[1]1304             hom(:,2,12,:) = SPREAD( zw, 2, statistic_regions+1 )
1305             IF ( prandtl_layer )  hom(nzb,2,12,:) = zu(1)
1306
1307          CASE ( 'w*u*' )
1308             dopr_index(i) = 13
[87]1309             dopr_unit(i)  = 'm2/s2'
[1]1310             hom(:,2,13,:) = SPREAD( zw, 2, statistic_regions+1 )
1311
1312          CASE ( 'w"v"' )
1313             dopr_index(i) = 14
[87]1314             dopr_unit(i)  = 'm2/s2'
[1]1315             hom(:,2,14,:) = SPREAD( zw, 2, statistic_regions+1 )
1316             IF ( prandtl_layer )  hom(nzb,2,14,:) = zu(1)
1317
1318          CASE ( 'w*v*' )
1319             dopr_index(i) = 15
[87]1320             dopr_unit(i)  = 'm2/s2'
[1]1321             hom(:,2,15,:) = SPREAD( zw, 2, statistic_regions+1 )
1322
1323          CASE ( 'w"pt"' )
1324             dopr_index(i) = 16
[87]1325             dopr_unit(i)  = 'K m/s'
[1]1326             hom(:,2,16,:) = SPREAD( zw, 2, statistic_regions+1 )
1327
1328          CASE ( 'w*pt*' )
1329             dopr_index(i) = 17
[87]1330             dopr_unit(i)  = 'K m/s'
[1]1331             hom(:,2,17,:) = SPREAD( zw, 2, statistic_regions+1 )
1332
1333          CASE ( 'wpt' )
1334             dopr_index(i) = 18
[87]1335             dopr_unit(i)  = 'K m/s'
[1]1336             hom(:,2,18,:) = SPREAD( zw, 2, statistic_regions+1 )
1337
1338          CASE ( 'wu' )
1339             dopr_index(i) = 19
[87]1340             dopr_unit(i)  = 'm2/s2'
[1]1341             hom(:,2,19,:) = SPREAD( zw, 2, statistic_regions+1 )
1342             IF ( prandtl_layer )  hom(nzb,2,19,:) = zu(1)
1343
1344          CASE ( 'wv' )
1345             dopr_index(i) = 20
[87]1346             dopr_unit(i)  = 'm2/s2'
[1]1347             hom(:,2,20,:) = SPREAD( zw, 2, statistic_regions+1 )
1348             IF ( prandtl_layer )  hom(nzb,2,20,:) = zu(1)
1349
1350          CASE ( 'w*pt*BC' )
1351             dopr_index(i) = 21
[87]1352             dopr_unit(i)  = 'K m/s'
[1]1353             hom(:,2,21,:) = SPREAD( zw, 2, statistic_regions+1 )
1354
1355          CASE ( 'wptBC' )
1356             dopr_index(i) = 22
[87]1357             dopr_unit(i)  = 'K m/s'
[1]1358             hom(:,2,22,:) = SPREAD( zw, 2, statistic_regions+1 )
1359
1360          CASE ( 'u*2' )
1361             dopr_index(i) = 30
[87]1362             dopr_unit(i)  = 'm2/s2'
[1]1363             hom(:,2,30,:) = SPREAD( zu, 2, statistic_regions+1 )
1364
1365          CASE ( 'v*2' )
1366             dopr_index(i) = 31
[87]1367             dopr_unit(i)  = 'm2/s2'
[1]1368             hom(:,2,31,:) = SPREAD( zu, 2, statistic_regions+1 )
1369
1370          CASE ( 'w*2' )
1371             dopr_index(i) = 32
[87]1372             dopr_unit(i)  = 'm2/s2'
[1]1373             hom(:,2,32,:) = SPREAD( zw, 2, statistic_regions+1 )
1374
1375          CASE ( 'pt*2' )
1376             dopr_index(i) = 33
[87]1377             dopr_unit(i)  = 'K2'
[1]1378             hom(:,2,33,:) = SPREAD( zu, 2, statistic_regions+1 )
1379
1380          CASE ( 'e*' )
1381             dopr_index(i) = 34
[87]1382             dopr_unit(i)  = 'm2/s2'
[1]1383             hom(:,2,34,:) = SPREAD( zu, 2, statistic_regions+1 )
1384
1385          CASE ( 'w*2pt*' )
1386             dopr_index(i) = 35
[87]1387             dopr_unit(i)  = 'K m2/s2'
[1]1388             hom(:,2,35,:) = SPREAD( zw, 2, statistic_regions+1 )
1389
1390          CASE ( 'w*pt*2' )
1391             dopr_index(i) = 36
[87]1392             dopr_unit(i)  = 'K2 m/s'
[1]1393             hom(:,2,36,:) = SPREAD( zw, 2, statistic_regions+1 )
1394
1395          CASE ( 'w*e*' )
1396             dopr_index(i) = 37
[87]1397             dopr_unit(i)  = 'm3/s3'
[1]1398             hom(:,2,37,:) = SPREAD( zw, 2, statistic_regions+1 )
1399
1400          CASE ( 'w*3' )
1401             dopr_index(i) = 38
[87]1402             dopr_unit(i)  = 'm3/s3'
[1]1403             hom(:,2,38,:) = SPREAD( zw, 2, statistic_regions+1 )
1404
1405          CASE ( 'Sw' )
1406             dopr_index(i) = 39
1407             hom(:,2,39,:) = SPREAD( zw, 2, statistic_regions+1 )
1408
1409          CASE ( 'q', '#q' )
1410             IF ( .NOT. cloud_physics )  THEN
1411                IF ( myid == 0 )  THEN
1412                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
1413                           data_output_pr(i),                          &
1414                           '    is not implemented for cloud_physics = FALSE'
1415                ENDIF
1416                CALL local_stop
1417             ELSE
1418                dopr_index(i) = 41
[87]1419                dopr_unit(i)  = 'kg/kg'
1420                hom(:,2,41,:) = SPREAD( zu, 2, statistic_regions+1 )
[1]1421                IF ( data_output_pr(i)(1:1) == '#' )  THEN
1422                   dopr_initial_index(i) = 26
1423                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
1424                   hom(nzb,2,26,:)       = 0.0    ! weil zu(nzb) negativ ist
1425                   data_output_pr(i)     = data_output_pr(i)(2:)
1426                ENDIF
1427             ENDIF
1428
1429          CASE ( 's', '#s' )
1430             IF ( .NOT. passive_scalar )  THEN
1431                IF ( myid == 0 )  THEN
1432                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
1433                           data_output_pr(i),                          &
1434                           '    is not implemented for passive_scalar = FALSE'
1435                ENDIF
1436                CALL local_stop
1437             ELSE
1438                dopr_index(i) = 41
[87]1439                dopr_unit(i)  = 'kg/m3'
1440                hom(:,2,41,:) = SPREAD( zu, 2, statistic_regions+1 )
[1]1441                IF ( data_output_pr(i)(1:1) == '#' )  THEN
1442                   dopr_initial_index(i) = 26
1443                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
1444                   hom(nzb,2,26,:)       = 0.0    ! weil zu(nzb) negativ ist
1445                   data_output_pr(i)     = data_output_pr(i)(2:)
1446                ENDIF
1447             ENDIF
1448
1449          CASE ( 'qv', '#qv' )
1450             IF ( .NOT. cloud_physics ) THEN
1451                dopr_index(i) = 41
[87]1452                dopr_unit(i)  = 'kg/kg'
1453                hom(:,2,41,:) = SPREAD( zu, 2, statistic_regions+1 )
[1]1454                IF ( data_output_pr(i)(1:1) == '#' )  THEN
1455                   dopr_initial_index(i) = 26
1456                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
1457                   hom(nzb,2,26,:)       = 0.0    ! weil zu(nzb) negativ ist
1458                   data_output_pr(i)     = data_output_pr(i)(2:)
1459                ENDIF
1460             ELSE
1461                dopr_index(i) = 42
[87]1462                dopr_unit(i)  = 'kg/kg'
1463                hom(:,2,42,:) = SPREAD( zu, 2, statistic_regions+1 )
[1]1464                IF ( data_output_pr(i)(1:1) == '#' )  THEN
1465                   dopr_initial_index(i) = 27
1466                   hom(:,2,27,:)         = SPREAD( zu, 2, statistic_regions+1 )
1467                   hom(nzb,2,27,:)       = 0.0    ! weil zu(nzb) negativ ist
1468                   data_output_pr(i)     = data_output_pr(i)(2:)
1469                ENDIF
1470             ENDIF
1471
1472          CASE ( 'lpt', '#lpt' )
1473             IF ( .NOT. cloud_physics ) THEN
1474                IF ( myid == 0 )  THEN
1475                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
1476                           data_output_pr(i),                          &
1477                           '    is not implemented for cloud_physics = FALSE'
1478                ENDIF
1479                CALL local_stop
1480             ELSE
1481                dopr_index(i) = 4
[87]1482                dopr_unit(i)  = 'K'
[1]1483                hom(:,2,4,:)  = SPREAD( zu, 2, statistic_regions+1 )
1484                IF ( data_output_pr(i)(1:1) == '#' )  THEN
1485                   dopr_initial_index(i) = 7
1486                   hom(:,2,7,:)          = SPREAD( zu, 2, statistic_regions+1 )
1487                   hom(nzb,2,7,:)        = 0.0    ! weil zu(nzb) negativ ist
1488                   data_output_pr(i)     = data_output_pr(i)(2:)
1489                ENDIF
1490             ENDIF
1491
1492          CASE ( 'vpt', '#vpt' )
1493             dopr_index(i) = 44
[87]1494             dopr_unit(i)  = 'K'
1495             hom(:,2,44,:) = SPREAD( zu, 2, statistic_regions+1 )
[1]1496             IF ( data_output_pr(i)(1:1) == '#' )  THEN
1497                dopr_initial_index(i) = 29
1498                hom(:,2,29,:)         = SPREAD( zu, 2, statistic_regions+1 )
1499                hom(nzb,2,29,:)       = 0.0    ! weil zu(nzb) negativ ist
1500                data_output_pr(i)     = data_output_pr(i)(2:)
1501             ENDIF
1502
1503          CASE ( 'w"vpt"' )
1504             dopr_index(i) = 45
[87]1505             dopr_unit(i)  = 'K m/s'
[1]1506             hom(:,2,45,:) = SPREAD( zw, 2, statistic_regions+1 )
1507
1508          CASE ( 'w*vpt*' )
1509             dopr_index(i) = 46
[87]1510             dopr_unit(i)  = 'K m/s'
[1]1511             hom(:,2,46,:) = SPREAD( zw, 2, statistic_regions+1 )
1512
1513          CASE ( 'wvpt' )
1514             dopr_index(i) = 47
[87]1515             dopr_unit(i)  = 'K m/s'
[1]1516             hom(:,2,47,:) = SPREAD( zw, 2, statistic_regions+1 )
1517
1518          CASE ( 'w"q"' )
1519             IF ( .NOT. cloud_physics ) THEN
1520                IF ( myid == 0 )  THEN
1521                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
1522                           data_output_pr(i),                          &
1523                           '    is not implemented for cloud_physics = FALSE'
1524                ENDIF
1525                CALL local_stop
1526             ELSE
1527                dopr_index(i) = 48
[87]1528                dopr_unit(i)  = 'kg/kg m/s'
[1]1529                hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
1530             ENDIF
1531
1532          CASE ( 'w*q*' )
1533             IF ( .NOT. cloud_physics ) THEN
1534                IF ( myid == 0 )  THEN
1535                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
1536                           data_output_pr(i),                          &
1537                           '    is not implemented for cloud_physics = FALSE'
1538                ENDIF
1539                CALL local_stop
1540             ELSE
1541                dopr_index(i) = 49
[87]1542                dopr_unit(i)  = 'kg/kg m/s'
[1]1543                hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
1544             ENDIF
1545
1546          CASE ( 'wq' )
1547             IF ( .NOT. cloud_physics ) THEN
1548                IF ( myid == 0 )  THEN
1549                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
1550                           data_output_pr(i),                          &
1551                           '    is not implemented for cloud_physics = FALSE'
1552                ENDIF
1553                CALL local_stop
1554             ELSE
1555                dopr_index(i) = 50
[87]1556                dopr_unit(i)  = 'kg/kg m/s'
[1]1557                hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
1558             ENDIF
1559
1560          CASE ( 'w"s"' )
1561             IF ( .NOT. passive_scalar ) THEN
1562                IF ( myid == 0 )  THEN
1563                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
1564                           data_output_pr(i),                          &
1565                           '    is not implemented for passive_scalar = FALSE'
1566                ENDIF
1567                CALL local_stop
1568             ELSE
1569                dopr_index(i) = 48
[87]1570                dopr_unit(i)  = 'kg/m3 m/s'
[1]1571                hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
1572             ENDIF
1573
1574          CASE ( 'w*s*' )
1575             IF ( .NOT. passive_scalar ) THEN
1576                IF ( myid == 0 )  THEN
1577                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
1578                           data_output_pr(i),                          &
1579                           '    is not implemented for passive_scalar = FALSE'
1580                ENDIF
1581                CALL local_stop
1582             ELSE
1583                dopr_index(i) = 49
[87]1584                dopr_unit(i)  = 'kg/m3 m/s'
[1]1585                hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
1586             ENDIF
1587
1588          CASE ( 'ws' )
1589             IF ( .NOT. passive_scalar ) THEN
1590                IF ( myid == 0 )  THEN
1591                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
1592                           data_output_pr(i),                          &
1593                           '    is not implemented for passive_scalar = FALSE'
1594                ENDIF
1595                CALL local_stop
1596             ELSE
1597                dopr_index(i) = 50
[87]1598                dopr_unit(i)  = 'kg/m3 m/s'
[1]1599                hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
1600             ENDIF
1601
1602          CASE ( 'w"qv"' )
[75]1603             IF ( humidity  .AND.  .NOT. cloud_physics ) &
[1]1604             THEN
1605                dopr_index(i) = 48
[87]1606                dopr_unit(i)  = 'kg/kg m/s'
[1]1607                hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
[75]1608             ELSEIF( humidity .AND. cloud_physics ) THEN
[1]1609                dopr_index(i) = 51
[87]1610                dopr_unit(i)  = 'kg/kg m/s'
[1]1611                hom(:,2,51,:) = SPREAD( zw, 2, statistic_regions+1 )
1612             ELSE
1613                IF ( myid == 0 )  THEN
1614                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
1615                           data_output_pr(i),                          &
1616                           '    is not implemented for cloud_physics = FALSE', &
[75]1617                           '    and                    humidity      = FALSE'
[1]1618                ENDIF
1619                CALL local_stop                   
1620             ENDIF
1621
1622          CASE ( 'w*qv*' )
[75]1623             IF ( humidity  .AND.  .NOT. cloud_physics ) &
[1]1624             THEN
1625                dopr_index(i) = 49
[87]1626                dopr_unit(i)  = 'kg/kg m/s'
[1]1627                hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
[75]1628             ELSEIF( humidity .AND. cloud_physics ) THEN
[1]1629                dopr_index(i) = 52
[87]1630                dopr_unit(i)  = 'kg/kg m/s'
[1]1631                hom(:,2,52,:) = SPREAD( zw, 2, statistic_regions+1 )
1632             ELSE
1633                IF ( myid == 0 )  THEN
1634                   PRINT*, '+++ check_parameters:  data_output_pr = ',         &
1635                           data_output_pr(i),                                  &
1636                           '    is not implemented for cloud_physics = FALSE', &
[75]1637                           '                       and humidity      = FALSE'
[1]1638                ENDIF
1639                CALL local_stop                   
1640             ENDIF
1641
1642          CASE ( 'wqv' )
[75]1643             IF ( humidity  .AND.  .NOT. cloud_physics ) &
[1]1644             THEN
1645                dopr_index(i) = 50
[87]1646                dopr_unit(i)  = 'kg/kg m/s'
[1]1647                hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
[75]1648             ELSEIF( humidity .AND. cloud_physics ) THEN
[1]1649                dopr_index(i) = 53
[87]1650                dopr_unit(i)  = 'kg/kg m/s'
[1]1651                hom(:,2,53,:) = SPREAD( zw, 2, statistic_regions+1 )
1652             ELSE
1653                IF ( myid == 0 )  THEN
1654                   PRINT*, '+++ check_parameters:  data_output_pr = ',         &
1655                           data_output_pr(i),                                  &
1656                           '    is not implemented for cloud_physics = FALSE', &
[75]1657                           '                       and humidity      = FALSE'
[1]1658                ENDIF
1659                CALL local_stop                   
1660             ENDIF
1661
1662          CASE ( 'ql' )
1663             IF ( .NOT. cloud_physics  .AND.  .NOT. cloud_droplets )  THEN
1664                IF ( myid == 0 )  THEN
1665                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
1666                           data_output_pr(i),                          &
1667                           '    is not implemented for cloud_physics = FALSE'
1668                ENDIF
1669                CALL local_stop
1670             ELSE
1671                dopr_index(i) = 54
[87]1672                dopr_unit(i)  = 'kg/kg'
[1]1673                hom(:,2,54,:)  = SPREAD( zu, 2, statistic_regions+1 )
1674             ENDIF
1675
1676          CASE ( 'w*u*u*/dz' )
1677             dopr_index(i) = 55
[87]1678             dopr_unit(i)  = 'm2/s3'
[1]1679             hom(:,2,55,:) = SPREAD( zu, 2, statistic_regions+1 )
1680
1681          CASE ( 'w*p*/dz' )
1682             dopr_index(i) = 56
[87]1683             dopr_unit(i)  = 'm2/s3'
[1]1684             hom(:,2,56,:) = SPREAD( zu, 2, statistic_regions+1 )
1685
1686          CASE ( 'w"e/dz' )
1687             dopr_index(i) = 57
[87]1688             dopr_unit(i)  = 'm2/s3'
[1]1689             hom(:,2,57,:) = SPREAD( zu, 2, statistic_regions+1 )
1690
1691          CASE ( 'u"pt"' )
1692             dopr_index(i) = 58
[87]1693             dopr_unit(i)  = 'K m/s'
[1]1694             hom(:,2,58,:) = SPREAD( zu, 2, statistic_regions+1 )
1695
1696          CASE ( 'u*pt*' )
1697             dopr_index(i) = 59
[87]1698             dopr_unit(i)  = 'K m/s'
[1]1699             hom(:,2,59,:) = SPREAD( zu, 2, statistic_regions+1 )
1700
1701          CASE ( 'upt_t' )
1702             dopr_index(i) = 60
[87]1703             dopr_unit(i)  = 'K m/s'
[1]1704             hom(:,2,60,:) = SPREAD( zu, 2, statistic_regions+1 )
1705
1706          CASE ( 'v"pt"' )
1707             dopr_index(i) = 61
[87]1708             dopr_unit(i)  = 'K m/s'
[1]1709             hom(:,2,61,:) = SPREAD( zu, 2, statistic_regions+1 )
1710             
1711          CASE ( 'v*pt*' )
1712             dopr_index(i) = 62
[87]1713             dopr_unit(i)  = 'K m/s'
[1]1714             hom(:,2,62,:) = SPREAD( zu, 2, statistic_regions+1 )
1715
1716          CASE ( 'vpt_t' )
1717             dopr_index(i) = 63
[87]1718             dopr_unit(i)  = 'K m/s'
[1]1719             hom(:,2,63,:) = SPREAD( zu, 2, statistic_regions+1 )
1720
1721
1722          CASE DEFAULT
[87]1723
1724             CALL user_check_data_output_pr( data_output_pr(i), i, unit )
1725
1726             IF ( unit == 'illegal' )  THEN
1727                IF ( myid == 0 )  THEN
1728                   IF ( data_output_pr_user(1) /= ' ' )  THEN
1729                      PRINT*, '+++ check_parameters:  illegal value for data_',&
1730                                   'output_pr or data_output_pr_user: "',      &
1731                                   TRIM( data_output_pr(i) ), '"'
1732                   ELSE
1733                      PRINT*, '+++ check_parameters:  illegal value for data_',&
1734                                   'output_pr: "', TRIM( data_output_pr(i) ),'"'
1735                   ENDIF
1736                ENDIF
1737                CALL local_stop
[1]1738             ENDIF
1739
1740       END SELECT
1741!
1742!--    Check to which of the predefined coordinate systems the profile belongs
1743       DO  k = 1, crmax
1744          IF ( INDEX( cross_profiles(k), ' '//TRIM( data_output_pr(i) )//' ' ) &
1745               /=0 ) &
1746          THEN
1747             dopr_crossindex(i) = k
1748             EXIT
1749          ENDIF
1750       ENDDO
1751!
1752!--    Generate the text for the labels of the PROFIL output file. "-characters
1753!--    must be substituted, otherwise PROFIL would interpret them as TeX
1754!--    control characters
1755       dopr_label(i) = data_output_pr(i)
1756       position = INDEX( dopr_label(i) , '"' )
1757       DO WHILE ( position /= 0 )
1758          dopr_label(i)(position:position) = ''''
1759          position = INDEX( dopr_label(i) , '"' )
1760       ENDDO
1761
1762    ENDDO
1763
1764!
1765!-- y-value range of the coordinate system (PROFIL).
1766!-- x-value range determined in plot_1d.
1767    cross_uymin = 0.0
1768    IF ( z_max_do1d == -1.0 )  THEN
1769       cross_uymax = zu(nzt+1)
1770    ELSEIF ( z_max_do1d < zu(nzb+1)  .OR.  z_max_do1d > zu(nzt+1) )  THEN
1771       IF ( myid == 0 )  PRINT*, '+++ check_parameters:  z_max_do1d=',  &
1772                                 z_max_do1d,' must be >= ', zu(nzb+1),  &
1773                                 ' or <= ', zu(nzt+1)
1774       CALL local_stop
1775    ELSE
1776       cross_uymax = z_max_do1d
1777    ENDIF
1778
1779!
1780!-- Check whether the chosen normalizing factor for the coordinate systems is
1781!-- permissible
1782    DO  i = 1, crmax
1783       SELECT CASE ( TRIM( cross_normalized_x(i) ) )  ! TRIM required on IBM
1784
1785          CASE ( '', 'wpt0', 'ws2', 'tsw2', 'ws3', 'ws2tsw', 'wstsw2' )
1786             j = 0
1787
1788          CASE DEFAULT
1789             IF ( myid == 0 )  THEN
1790                PRINT*, '+++ check_parameters: unknown normalize method'
1791                PRINT*, '    cross_normalized_x="',cross_normalized_x(i),'"'
1792             ENDIF
1793             CALL local_stop
1794
1795       END SELECT
1796       SELECT CASE ( TRIM( cross_normalized_y(i) ) )  ! TRIM required on IBM
1797
1798          CASE ( '', 'z_i' )
1799             j = 0
1800
1801          CASE DEFAULT
1802             IF ( myid == 0 )  THEN
1803                PRINT*, '+++ check_parameters: unknown normalize method'
1804                PRINT*, '    cross_normalized_y="',cross_normalized_y(i),'"'
1805             ENDIF
1806             CALL local_stop
1807
1808       END SELECT
1809    ENDDO
1810!
1811!-- Check normalized y-value range of the coordinate system (PROFIL)
1812    IF ( z_max_do1d_normalized /= -1.0  .AND.  z_max_do1d_normalized <= 0.0 ) &
1813    THEN
1814       IF ( myid == 0 )  PRINT*,'+++ check_parameters:  z_max_do1d_normalize', &
1815                                'd=', z_max_do1d_normalized, ' must be >= 0.0'
1816       CALL local_stop
1817    ENDIF
1818
1819
1820!
1821!-- Append user-defined data output variables to the standard data output
1822    IF ( data_output_user(1) /= ' ' )  THEN
1823       i = 1
1824       DO  WHILE ( data_output(i) /= ' '  .AND.  i <= 100 )
1825          i = i + 1
1826       ENDDO
1827       j = 1
1828       DO  WHILE ( data_output_user(j) /= ' '  .AND.  j <= 100 )
1829          IF ( i > 100 )  THEN
1830             IF ( myid == 0 )  THEN
1831                PRINT*, '+++ check_parameters: number of output quantitities', &
1832                             ' given by data_output and data_output_user'
1833                PRINT*, '                      exceeds the limit of 100'
1834             ENDIF
1835             CALL local_stop
1836          ENDIF
1837          data_output(i) = data_output_user(j)
1838          i = i + 1
1839          j = j + 1
1840       ENDDO
1841    ENDIF
1842
1843!
1844!-- Check and set steering parameters for 2d/3d data output and averaging
1845    i   = 1
1846    DO  WHILE ( data_output(i) /= ' '  .AND.  i <= 100 )
1847!
1848!--    Check for data averaging
1849       ilen = LEN_TRIM( data_output(i) )
1850       j = 0                                                 ! no data averaging
1851       IF ( ilen > 3 )  THEN
1852          IF ( data_output(i)(ilen-2:ilen) == '_av' )  THEN
1853             j = 1                                           ! data averaging
1854             data_output(i) = data_output(i)(1:ilen-3)
1855          ENDIF
1856       ENDIF
1857!
1858!--    Check for cross section or volume data
1859       ilen = LEN_TRIM( data_output(i) )
1860       k = 0                                                   ! 3d data
1861       var = data_output(i)(1:ilen)
1862       IF ( ilen > 3 )  THEN
1863          IF ( data_output(i)(ilen-2:ilen) == '_xy'  .OR. &
1864               data_output(i)(ilen-2:ilen) == '_xz'  .OR. &
1865               data_output(i)(ilen-2:ilen) == '_yz' )  THEN
1866             k = 1                                             ! 2d data
1867             var = data_output(i)(1:ilen-3)
1868          ENDIF
1869       ENDIF
1870!
1871!--    Check for allowed value and set units
1872       SELECT CASE ( TRIM( var ) )
1873
1874          CASE ( 'e' )
1875             IF ( constant_diffusion )  THEN
1876                IF ( myid == 0 )  THEN
1877                   PRINT*, '+++ check_parameters: output of "', TRIM( var ), &
1878                                '" requires constant_diffusion = .FALSE.'
1879                ENDIF
1880                CALL local_stop
1881             ENDIF
1882             unit = 'm2/s2'
1883
1884          CASE ( 'pc', 'pr' )
1885             IF ( .NOT. particle_advection )  THEN
1886                IF ( myid == 0 )  THEN
1887                   PRINT*, '+++ check_parameters: output of "', TRIM( var ), &
1888                                '" requires particle package'
1889                   PRINT*, '                      (mrun-option "-p particles")'
1890                ENDIF
1891                CALL local_stop
1892             ENDIF
1893             IF ( TRIM( var ) == 'pc' )  unit = 'number'
1894             IF ( TRIM( var ) == 'pr' )  unit = 'm'
1895
1896          CASE ( 'q', 'vpt' )
[75]1897             IF ( .NOT. humidity )  THEN
[1]1898                IF ( myid == 0 )  THEN
1899                   PRINT*, '+++ check_parameters: output of "', TRIM( var ), &
[75]1900                                '" requires humidity = .TRUE.'
[1]1901                ENDIF
1902                CALL local_stop
1903             ENDIF
1904             IF ( TRIM( var ) == 'q'   )  unit = 'kg/kg'
1905             IF ( TRIM( var ) == 'vpt' )  unit = 'K'
1906
1907          CASE ( 'ql' )
1908             IF ( .NOT. ( cloud_physics  .OR.  cloud_droplets ) )  THEN
1909                IF ( myid == 0 )  THEN
1910                   PRINT*, '+++ check_parameters: output of "', TRIM( var ), &
1911                                '" requires cloud_physics = .TRUE.'
1912                   PRINT*, '                      or cloud_droplets = .TRUE.'
1913                ENDIF
1914                CALL local_stop
1915             ENDIF
1916             unit = 'kg/kg'
1917
1918          CASE ( 'ql_c', 'ql_v', 'ql_vp' )
1919             IF ( .NOT. cloud_droplets )  THEN
1920                IF ( myid == 0 )  THEN
1921                   PRINT*, '+++ check_parameters: output of "', TRIM( var ), &
1922                                '" requires cloud_droplets = .TRUE.'
1923                ENDIF
1924                CALL local_stop
1925             ENDIF
1926             IF ( TRIM( var ) == 'ql_c'  )  unit = 'kg/kg'
1927             IF ( TRIM( var ) == 'ql_v'  )  unit = 'm3'
1928             IF ( TRIM( var ) == 'ql_vp' )  unit = 'none'
1929
1930          CASE ( 'qv' )
1931             IF ( .NOT. cloud_physics )  THEN
1932                IF ( myid == 0 )  THEN
1933                   PRINT*, '+++ check_parameters: output of "', TRIM( var ), &
1934                                '" requires cloud_physics = .TRUE.'
1935                ENDIF
1936                CALL local_stop
1937             ENDIF
1938             unit = 'kg/kg'
1939
1940          CASE ( 's' )
1941             IF ( .NOT. passive_scalar )  THEN
1942                IF ( myid == 0 )  THEN
1943                   PRINT*, '+++ check_parameters: output of "', TRIM( var ), &
1944                                '" requires passive_scalar = .TRUE.'
1945                ENDIF
1946                CALL local_stop
1947             ENDIF
1948             unit = 'conc'
1949
[72]1950          CASE ( 'u*', 't*', 'lwp*', 'pra*', 'prr*', 'z0*' )
[1]1951             IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
1952                IF ( myid == 0 )  THEN
1953                   PRINT*, '+++ check_parameters:  illegal value for data_',&
1954                                'output: "', TRIM( var ), '" is only allowed'
1955                   PRINT*, '                       for horizontal cross section'
1956                ENDIF
1957                CALL local_stop
1958             ENDIF
1959             IF ( TRIM( var ) == 'lwp*'  .AND.  .NOT. cloud_physics )  THEN
1960                IF ( myid == 0 )  THEN
1961                   PRINT*, '+++ check_parameters: output of "', TRIM( var ), &
1962                                '" requires cloud_physics = .TRUE.'
1963                ENDIF
1964                CALL local_stop
1965             ENDIF
[72]1966             IF ( TRIM( var ) == 'pra*'  .AND.  .NOT. precipitation )  THEN
1967                IF ( myid == 0 )  THEN
1968                   PRINT*, '+++ check_parameters: output of "', TRIM( var ), &
1969                                '" requires precipitation = .TRUE.'
1970                ENDIF
1971                CALL local_stop
1972             ENDIF
1973             IF ( TRIM( var ) == 'pra*'  .AND.  j == 1 )  THEN
1974                IF ( myid == 0 )  THEN
1975                   PRINT*, '+++ check_parameters: temporal averaging of ', &
1976                           ' precipitation amount "', TRIM( var ),         &
1977                           '" not possible' 
1978                ENDIF
1979                CALL local_stop
1980             ENDIF
1981             IF ( TRIM( var ) == 'prr*'  .AND.  .NOT. precipitation )  THEN
1982                IF ( myid == 0 )  THEN
1983                   PRINT*, '+++ check_parameters: output of "', TRIM( var ), &
1984                                '" requires precipitation = .TRUE.'
1985                ENDIF
1986                CALL local_stop
1987             ENDIF
1988
1989
[1]1990             IF ( TRIM( var ) == 'u*'   )  unit = 'm/s'
1991             IF ( TRIM( var ) == 't*'   )  unit = 'K'
1992             IF ( TRIM( var ) == 'lwp*' )  unit = 'kg/kg*m'
[72]1993             IF ( TRIM( var ) == 'pra*' )  unit = 'mm'
1994             IF ( TRIM( var ) == 'prr*' )  unit = 'mm/s'
1995             IF ( TRIM( var ) == 'z0*'  )  unit = 'm'
[1]1996
1997          CASE ( 'p', 'pt', 'u', 'v', 'w' )
1998             IF ( TRIM( var ) == 'p'  )  unit = 'Pa'
1999             IF ( TRIM( var ) == 'pt' )  unit = 'K'
2000             IF ( TRIM( var ) == 'u'  )  unit = 'm/s'
2001             IF ( TRIM( var ) == 'v'  )  unit = 'm/s'
2002             IF ( TRIM( var ) == 'w'  )  unit = 'm/s'
2003             CONTINUE
2004
2005          CASE DEFAULT
2006             CALL user_check_data_output( var, unit )
2007
2008             IF ( unit == 'illegal' )  THEN
2009                IF ( myid == 0 )  THEN
2010                   IF ( data_output_user(1) /= ' ' )  THEN
2011                      PRINT*, '+++ check_parameters:  illegal value for data_',&
2012                                   'output or data_output_user: "',            &
2013                                   TRIM( data_output(i) ), '"'
2014                   ELSE
2015                      PRINT*, '+++ check_parameters:  illegal value for data_',&
2016                                   'output: "', TRIM( data_output(i) ), '"'
2017                   ENDIF
2018                ENDIF
2019                CALL local_stop
2020             ENDIF
2021
2022       END SELECT
2023!
2024!--    Set the internal steering parameters appropriately
2025       IF ( k == 0 )  THEN
2026          do3d_no(j)              = do3d_no(j) + 1
2027          do3d(j,do3d_no(j))      = data_output(i)
2028          do3d_unit(j,do3d_no(j)) = unit
2029       ELSE
2030          do2d_no(j)              = do2d_no(j) + 1
2031          do2d(j,do2d_no(j))      = data_output(i)
2032          do2d_unit(j,do2d_no(j)) = unit
2033          IF ( data_output(i)(ilen-2:ilen) == '_xy' )  THEN
2034             data_output_xy(j) = .TRUE.
2035          ENDIF
2036          IF ( data_output(i)(ilen-2:ilen) == '_xz' )  THEN
2037             data_output_xz(j) = .TRUE.
2038          ENDIF
2039          IF ( data_output(i)(ilen-2:ilen) == '_yz' )  THEN
2040             data_output_yz(j) = .TRUE.
2041          ENDIF
2042       ENDIF
2043
2044       IF ( j == 1 )  THEN
2045!
2046!--       Check, if variable is already subject to averaging
2047          found = .FALSE.
2048          DO  k = 1, doav_n
2049             IF ( TRIM( doav(k) ) == TRIM( var ) )  found = .TRUE.
2050          ENDDO
2051
2052          IF ( .NOT. found )  THEN
2053             doav_n = doav_n + 1
2054             doav(doav_n) = var
2055          ENDIF
2056       ENDIF
2057
2058       i = i + 1
2059    ENDDO
2060
2061!
2062!-- Store sectional planes in one shared array
2063    section(:,1) = section_xy
2064    section(:,2) = section_xz
2065    section(:,3) = section_yz
2066
2067!
2068!-- Upper plot limit (grid point value) for 1D profiles
2069    IF ( z_max_do1d == -1.0 )  THEN
2070       nz_do1d = nzt+1
2071    ELSE
2072       DO  k = nzb+1, nzt+1
2073          nz_do1d = k
2074          IF ( zw(k) > z_max_do1d )  EXIT
2075       ENDDO
2076    ENDIF
2077
2078!
2079!-- Upper plot limit for 2D vertical sections
2080    IF ( z_max_do2d == -1.0 )  z_max_do2d = zu(nzt)
2081    IF ( z_max_do2d < zu(nzb+1)  .OR.  z_max_do2d > zu(nzt) )  THEN
2082       IF ( myid == 0 )  PRINT*, '+++ check_parameters:  z_max_do2d=',        &
2083                                 z_max_do2d, ' must be >= ', zu(nzb+1),       &
2084                                 '(zu(nzb+1)) and <= ', zu(nzt), ' (zu(nzt))'
2085       CALL local_stop
2086    ENDIF
2087
2088!
2089!-- Upper plot limit for 3D arrays
2090    IF ( nz_do3d == -9999 )  nz_do3d = nzt + 1
2091
2092!
2093!-- Determine and check accuracy for compressed 3D plot output
2094    IF ( do3d_compress )  THEN
2095!
2096!--    Compression only permissible on T3E machines
2097       IF ( host(1:3) /= 't3e' )  THEN
2098          IF ( myid == 0 )  THEN
2099             PRINT*, '+++ check_parameters: do3d_compress = .TRUE. not allow', &
2100                          'ed on host "', TRIM( host ), '"'
2101          ENDIF
2102          CALL local_stop
2103       ENDIF
2104
2105       i = 1
2106       DO  WHILE ( do3d_comp_prec(i) /= ' ' )
2107
2108          ilen = LEN_TRIM( do3d_comp_prec(i) )
2109          IF ( LLT( do3d_comp_prec(i)(ilen:ilen), '0' ) .OR. &
2110               LGT( do3d_comp_prec(i)(ilen:ilen), '9' ) )  THEN
2111             IF ( myid == 0 )  THEN
2112                PRINT*, '+++ check_parameters: illegal precision: ', &
2113                        'do3d_comp_prec(', i, ')="', TRIM(do3d_comp_prec(i)),'"'
2114             ENDIF
2115             CALL local_stop
2116          ENDIF
2117
2118          prec = IACHAR( do3d_comp_prec(i)(ilen:ilen) ) - IACHAR( '0' )
2119          var = do3d_comp_prec(i)(1:ilen-1)
2120
2121          SELECT CASE ( var )
2122
2123             CASE ( 'u' )
2124                j = 1
2125             CASE ( 'v' )
2126                j = 2
2127             CASE ( 'w' )
2128                j = 3
2129             CASE ( 'p' )
2130                j = 4
2131             CASE ( 'pt' )
2132                j = 5
2133
2134             CASE DEFAULT
2135                IF ( myid == 0 )  THEN
2136                   PRINT*, '+++ check_parameters: unknown variable in ', &
2137                               'assignment'
2138                   PRINT*, '    do3d_comp_prec(', i, ')="', &
2139                           TRIM( do3d_comp_prec(i) ),'"'
2140                ENDIF
2141                CALL local_stop               
2142
2143          END SELECT
2144
2145          plot_3d_precision(j)%precision = prec
2146          i = i + 1
2147
2148       ENDDO
2149    ENDIF
2150
2151!
2152!-- Check the data output format(s)
2153    IF ( data_output_format(1) == ' ' )  THEN
2154!
2155!--    Default value
2156       netcdf_output = .TRUE.
2157    ELSE
2158       i = 1
2159       DO  WHILE ( data_output_format(i) /= ' ' )
2160
2161          SELECT CASE ( data_output_format(i) )
2162
2163             CASE ( 'netcdf' )
2164                netcdf_output = .TRUE.
2165             CASE ( 'iso2d' )
2166                iso2d_output  = .TRUE.
2167             CASE ( 'profil' )
2168                profil_output = .TRUE.
2169             CASE ( 'avs' )
2170                avs_output    = .TRUE.
2171
2172             CASE DEFAULT
2173                IF ( myid == 0 )  THEN
2174                   PRINT*, '+++ check_parameters:'
2175                   PRINT*, '    unknown value for data_output_format "', &
2176                                TRIM( data_output_format(i) ),'"'
2177                ENDIF
2178                CALL local_stop               
2179
2180          END SELECT
2181
2182          i = i + 1
2183          IF ( i > 10 )  EXIT
2184
2185       ENDDO
2186
2187    ENDIF
2188
2189!
2190!-- Check netcdf precison
2191    ldum = .FALSE.
2192    CALL define_netcdf_header( 'ch', ldum, 0 )
2193
2194!
2195!-- Check, whether a constant diffusion coefficient shall be used
2196    IF ( km_constant /= -1.0 )  THEN
2197       IF ( km_constant < 0.0 )  THEN
2198          IF ( myid == 0 )  PRINT*, '+++ check_parameters:  km_constant=', &
2199                                    km_constant, ' < 0.0'
2200          CALL local_stop
2201       ELSE
2202          IF ( prandtl_number < 0.0 )  THEN
2203             IF ( myid == 0 )  PRINT*,'+++ check_parameters:  prandtl_number=',&
2204                                      prandtl_number, ' < 0.0'
2205             CALL local_stop
2206          ENDIF
2207          constant_diffusion = .TRUE.
2208
2209          IF ( prandtl_layer )  THEN
2210             IF ( myid == 0 )  PRINT*, '+++ check_parameters:  prandtl_layer ',&
2211                                       'is not allowed with fixed value of km'
2212             CALL local_stop
2213          ENDIF
2214       ENDIF
2215    ENDIF
2216
2217!
2218!-- In case of non-cyclic lateral boundaries, set the default maximum value
2219!-- for the horizontal diffusivity used within the outflow damping layer,
2220!-- and check/set the width of the damping layer
2221    IF ( bc_lr /= 'cyclic' ) THEN
2222       IF ( km_damp_max == -1.0 )  THEN
2223          km_damp_max = 0.5 * dx
2224       ENDIF
2225       IF ( outflow_damping_width == -1.0 )  THEN
2226          outflow_damping_width = MIN( 20, nx/2 )
2227       ENDIF
2228       IF ( outflow_damping_width <= 0  .OR.  outflow_damping_width > nx )  THEN
2229          IF ( myid == 0 )  PRINT*, '+++ check_parameters: outflow_damping w',&
2230                                    'idth out of range'
2231          CALL local_stop
2232       ENDIF
2233    ENDIF
2234
2235    IF ( bc_ns /= 'cyclic' )  THEN
2236       IF ( km_damp_max == -1.0 )  THEN
2237          km_damp_max = 0.5 * dy
2238       ENDIF
2239       IF ( outflow_damping_width == -1.0 )  THEN
2240          outflow_damping_width = MIN( 20, ny/2 )
2241       ENDIF
2242       IF ( outflow_damping_width <= 0  .OR.  outflow_damping_width > ny )  THEN
2243          IF ( myid == 0 )  PRINT*, '+++ check_parameters: outflow_damping w',&
2244                                    'idth out of range'
2245          CALL local_stop
2246       ENDIF
2247    ENDIF
2248
2249!
2250!-- Check value range for rif
2251    IF ( rif_min >= rif_max )  THEN
2252       IF ( myid == 0 )  PRINT*, '+++ check_parameters:  rif_min=', rif_min, &
2253                                 ' must be less than rif_max=', rif_max
2254       CALL local_stop
2255    ENDIF
2256
2257!
2258!-- Determine upper and lower hight level indices for random perturbations
2259    IF ( disturbance_level_b == -1.0 )  THEN
2260       disturbance_level_b     = zu(nzb+3)
2261       disturbance_level_ind_b = nzb + 3
2262    ELSEIF ( disturbance_level_b < zu(3) )  THEN
2263       IF ( myid == 0 )  PRINT*, '+++ check_parameters:  disturbance_level_b=',&
2264                                 disturbance_level_b, ' must be >= ',zu(3),    &
2265                                 '(zu(3))'
2266       CALL local_stop
2267    ELSEIF ( disturbance_level_b > zu(nzt-2) )  THEN
2268       IF ( myid == 0 )  PRINT*, '+++ check_parameters:  disturbance_level_b=',&
2269                                 disturbance_level_b, ' must be <= ',zu(nzt-2),&
2270                                 '(zu(nzt-2))'
2271       CALL local_stop
2272    ELSE
2273       DO  k = 3, nzt-2
2274          IF ( disturbance_level_b <= zu(k) )  THEN
2275             disturbance_level_ind_b = k
2276             EXIT
2277          ENDIF
2278       ENDDO
2279    ENDIF
2280
2281    IF ( disturbance_level_t == -1.0 )  THEN
2282       disturbance_level_t     = zu(nzt/3)
2283       disturbance_level_ind_t = nzt / 3
2284    ELSEIF ( disturbance_level_t > zu(nzt-2) )  THEN
2285       IF ( myid == 0 )  PRINT*, '+++ check_parameters:  disturbance_level_t=',&
2286                                 disturbance_level_t, ' must be <= ',zu(nzt-2),&
2287                                 '(zu(nzt-2))'
2288       CALL local_stop
2289    ELSEIF ( disturbance_level_t < disturbance_level_b )  THEN
2290       IF ( myid == 0 )  PRINT*, '+++ check_parameters:  disturbance_level_t=',&
2291                                 disturbance_level_t, ' must be >= ',          &
2292                                 'disturbance_level_b=', disturbance_level_b
2293       CALL local_stop
2294    ELSE
2295       DO  k = 3, nzt-2
2296          IF ( disturbance_level_t <= zu(k) )  THEN
2297             disturbance_level_ind_t = k
2298             EXIT
2299          ENDIF
2300       ENDDO
2301    ENDIF
2302
2303!
2304!-- Check again whether the levels determined this way are ok.
2305!-- Error may occur at automatic determination and too few grid points in
2306!-- z-direction.
2307    IF ( disturbance_level_ind_t < disturbance_level_ind_b )  THEN
2308       IF ( myid == 0 )  PRINT*, '+++ check_parameters:  ',                &
2309                                 'disturbance_level_ind_t=',               &
2310                                 disturbance_level_ind_t, ' must be >= ',  &
2311                                 'disturbance_level_ind_b=',               &
2312                                 disturbance_level_ind_b
2313       CALL local_stop
2314    ENDIF
2315
2316!
2317!-- Determine the horizontal index range for random perturbations.
2318!-- In case of non-cyclic horizontal boundaries, no perturbations are imposed
2319!-- near the inflow and the perturbation area is further limited to ...(1)
2320!-- after the initial phase of the flow.
2321    dist_nxl = 0;  dist_nxr = nx
2322    dist_nys = 0;  dist_nyn = ny
2323    IF ( bc_lr /= 'cyclic' )  THEN
2324       IF ( inflow_disturbance_begin == -1 )  THEN
2325          inflow_disturbance_begin = MIN( 10, nx/2 )
2326       ENDIF
2327       IF ( inflow_disturbance_begin < 0  .OR.  inflow_disturbance_begin > nx )&
2328       THEN
2329          IF ( myid == 0 )  PRINT*, '+++ check_parameters: inflow_disturbance',&
2330                                    '_begin out of range'
2331          CALL local_stop
2332       ENDIF
2333       IF ( inflow_disturbance_end == -1 )  THEN
2334          inflow_disturbance_end = MIN( 100, 3*nx/4 )
2335       ENDIF
2336       IF ( inflow_disturbance_end < 0  .OR.  inflow_disturbance_end > nx )    &
2337       THEN
2338          IF ( myid == 0 )  PRINT*, '+++ check_parameters: inflow_disturbance',&
2339                                    '_end out of range'
2340          CALL local_stop
2341       ENDIF
2342    ELSEIF ( bc_ns /= 'cyclic' )  THEN
2343       IF ( inflow_disturbance_begin == -1 )  THEN
2344          inflow_disturbance_begin = MIN( 10, ny/2 )
2345       ENDIF
2346       IF ( inflow_disturbance_begin < 0  .OR.  inflow_disturbance_begin > ny )&
2347       THEN
2348          IF ( myid == 0 )  PRINT*, '+++ check_parameters: inflow_disturbance',&
2349                                    '_begin out of range'
2350          CALL local_stop
2351       ENDIF
2352       IF ( inflow_disturbance_end == -1 )  THEN
2353          inflow_disturbance_end = MIN( 100, 3*ny/4 )
2354       ENDIF
2355       IF ( inflow_disturbance_end < 0  .OR.  inflow_disturbance_end > ny )    &
2356       THEN
2357          IF ( myid == 0 )  PRINT*, '+++ check_parameters: inflow_disturbance',&
2358                                    '_end out of range'
2359          CALL local_stop
2360       ENDIF
2361    ENDIF
2362
[73]2363    IF ( bc_lr == 'radiation/dirichlet' )  THEN
[1]2364       dist_nxr    = nx - inflow_disturbance_begin
2365       dist_nxl(1) = nx - inflow_disturbance_end
[73]2366    ELSEIF ( bc_lr == 'dirichlet/radiation' )  THEN
[1]2367       dist_nxl    = inflow_disturbance_begin
2368       dist_nxr(1) = inflow_disturbance_end
[73]2369    ENDIF
2370    IF ( bc_ns == 'dirichlet/radiation' )  THEN
[1]2371       dist_nyn    = ny - inflow_disturbance_begin
2372       dist_nys(1) = ny - inflow_disturbance_end
[73]2373    ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
[1]2374       dist_nys    = inflow_disturbance_begin
2375       dist_nyn(1) = inflow_disturbance_end
2376    ENDIF
2377
2378!
2379!-- Check random generator
2380    IF ( random_generator /= 'system-specific'  .AND. &
2381         random_generator /= 'numerical-recipes' )  THEN
2382       IF ( myid == 0 )  THEN
2383          PRINT*, '+++ check_parameters:'
2384          PRINT*, '    unknown random generator: random_generator=', &
2385                  random_generator
2386       ENDIF
2387       CALL local_stop
2388    ENDIF
2389
2390!
2391!-- Determine damping level index for 1D model
2392    IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
2393       IF ( damp_level_1d == -1.0 )  THEN
2394          damp_level_1d     = zu(nzt+1)
2395          damp_level_ind_1d = nzt + 1
2396       ELSEIF ( damp_level_1d < 0.0  .OR.  damp_level_1d > zu(nzt+1) )  THEN
2397          IF ( myid == 0 )  PRINT*, '+++ check_parameters:  damp_level_1d=', &
2398                                    damp_level_1d, ' must be > 0.0 and < ',  &
2399                                    'zu(nzt+1)', zu(nzt+1)
2400          CALL local_stop
2401       ELSE
2402          DO  k = 1, nzt+1
2403             IF ( damp_level_1d <= zu(k) )  THEN
2404                damp_level_ind_1d = k
2405                EXIT
2406             ENDIF
2407          ENDDO
2408       ENDIF
2409    ENDIF
2410!
2411!-- Check some other 1d-model parameters
2412    IF ( TRIM( mixing_length_1d ) /= 'as_in_3d_model'  .AND. &
2413         TRIM( mixing_length_1d ) /= 'blackadar' )  THEN
2414       IF ( myid == 0 )  PRINT*, '+++ check_parameters: mixing_length_1d = "', &
2415                                 TRIM( mixing_length_1d ), '" is unknown'
2416       CALL local_stop
2417    ENDIF
2418    IF ( TRIM( dissipation_1d ) /= 'as_in_3d_model'  .AND. &
2419         TRIM( dissipation_1d ) /= 'detering' )  THEN
2420       IF ( myid == 0 )  PRINT*, '+++ check_parameters: dissipation_1d = "', &
2421                                 TRIM( dissipation_1d ), '" is unknown'
2422       CALL local_stop
2423    ENDIF
2424
2425!
2426!-- Set time for the next user defined restart (time_restart is the
2427!-- internal parameter for steering restart events)
2428    IF ( restart_time /= 9999999.9 )  THEN
2429       IF ( restart_time > simulated_time )  time_restart = restart_time
2430    ELSE
2431!
2432!--    In case of a restart run, set internal parameter to default (no restart)
2433!--    if the NAMELIST-parameter restart_time is omitted
2434       time_restart = 9999999.9
2435    ENDIF
2436
2437!
2438!-- Set default value of the time needed to terminate a model run
2439    IF ( termination_time_needed == -1.0 )  THEN
2440       IF ( host(1:3) == 'ibm' )  THEN
2441          termination_time_needed = 300.0
2442       ELSE
2443          termination_time_needed = 35.0
2444       ENDIF
2445    ENDIF
2446
2447!
2448!-- Check the time needed to terminate a model run
2449    IF ( host(1:3) == 't3e' )  THEN
2450!
2451!--    Time needed must be at least 30 seconds on all CRAY machines, because
2452!--    MPP_TREMAIN gives the remaining CPU time only in steps of 30 seconds
2453       IF ( termination_time_needed <= 30.0 )  THEN
2454          IF ( myid == 0 )  THEN
2455             PRINT*, '+++ check_parameters:  termination_time_needed', &
2456                      termination_time_needed
2457             PRINT*, '                       must be > 30.0 on host "', host, &
2458                     '"'
2459          ENDIF
2460          CALL local_stop
2461       ENDIF
2462    ELSEIF ( host(1:3) == 'ibm' )  THEN
2463!
2464!--    On IBM-regatta machines the time should be at least 300 seconds,
2465!--    because the job time consumed before executing palm (for compiling,
2466!--    copying of files, etc.) has to be regarded
2467       IF ( termination_time_needed < 300.0 )  THEN
2468          IF ( myid == 0 )  THEN
2469             PRINT*, '+++ WARNING: check_parameters:  termination_time_',  &
2470                         'needed', termination_time_needed
2471             PRINT*, '                                should be >= 300.0', &
2472                         ' on host "', host, '"'
2473          ENDIF
2474       ENDIF
2475    ENDIF
2476
2477
2478 END SUBROUTINE check_parameters
Note: See TracBrowser for help on using the repository browser.