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

Last change on this file since 177 was 177, checked in by steinfeld, 16 years ago

Bug in the construction of initial vertical profiles fixed

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