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

Last change on this file since 182 was 181, checked in by raasch, 16 years ago

bugfixes + adjustments for SGI ICE system

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