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

Last change on this file since 1455 was 1455, checked in by heinze, 10 years ago

empty time records in volume, cross-section and masked data output prevented in case of non-parallel netcdf-output in restart runs

  • Property svn:keywords set to Id
File size: 151.0 KB
Line 
1 SUBROUTINE check_parameters
2
3!--------------------------------------------------------------------------------!
4! This file is part of PALM.
5!
6! PALM is free software: you can redistribute it and/or modify it under the terms
7! of the GNU General Public License as published by the Free Software Foundation,
8! either version 3 of the License, or (at your option) any later version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2014 Leibniz Universitaet Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22! empty time records in volume, cross-section and masked data output prevented 
23! in case of non-parallel netcdf-output in restart runs
24!
25! Former revisions:
26! -----------------
27! $Id: check_parameters.f90 1455 2014-08-29 10:47:47Z heinze $
28!
29! 1429 2014-07-15 12:53:45Z knoop
30! run_description_header exended to provide ensemble_member_nr if specified
31!
32! 1425 2014-07-05 10:57:53Z knoop
33! bugfix: perturbation domain modified for parallel random number generator
34!
35! 1402 2014-05-09 14:25:13Z raasch
36! location messages modified
37!
38! 1400 2014-05-09 14:03:54Z knoop
39! Check random generator extended by option random-parallel
40!
41! 1384 2014-05-02 14:31:06Z raasch
42! location messages added
43!
44! 1365 2014-04-22 15:03:56Z boeske
45! Usage of large scale forcing for pt and q enabled:
46! output for profiles of large scale advection (td_lsa_lpt, td_lsa_q),
47! large scale subsidence (td_sub_lpt, td_sub_q)
48! and nudging tendencies (td_nud_lpt, td_nud_q, td_nud_u and td_nud_v) added,
49! check if use_subsidence_tendencies is used correctly
50!
51! 1361 2014-04-16 15:17:48Z hoffmann
52! PA0363 removed
53! PA0362 changed
54!
55! 1359 2014-04-11 17:15:14Z hoffmann
56! Do not allow the execution of PALM with use_particle_tails, since particle
57! tails are currently not supported by our new particle structure.
58!
59! PA0084 not necessary for new particle structure
60!
61! 1353 2014-04-08 15:21:23Z heinze
62! REAL constants provided with KIND-attribute
63!
64! 1330 2014-03-24 17:29:32Z suehring
65! In case of SGS-particle velocity advection of TKE is also allowed with
66! dissipative 5th-order scheme.
67!
68! 1327 2014-03-21 11:00:16Z raasch
69! "baroclinicity" renamed "baroclinity", "ocean version" replaced by "ocean mode"
70! bugfix: duplicate error message 56 removed,
71! check of data_output_format and do3d_compress removed
72!
73! 1322 2014-03-20 16:38:49Z raasch
74! some REAL constants defined as wp-kind
75!
76! 1320 2014-03-20 08:40:49Z raasch
77! Kind-parameters added to all INTEGER and REAL declaration statements,
78! kinds are defined in new module kinds,
79! revision history before 2012 removed,
80! comment fields (!:) to be used for variable explanations added to
81! all variable declaration statements
82!
83! 1308 2014-03-13 14:58:42Z fricke
84! +netcdf_data_format_save
85! Calculate fixed number of output time levels for parallel netcdf output.
86! For masked data, parallel netcdf output is not tested so far, hence
87! netcdf_data_format is switched back to non-paralell output.
88!
89! 1299 2014-03-06 13:15:21Z heinze
90! enable usage of large_scale subsidence in combination with large_scale_forcing
91! output for profile of large scale vertical velocity w_subs added
92!
93! 1276 2014-01-15 13:40:41Z heinze
94! Use LSF_DATA also in case of Dirichlet bottom boundary condition for scalars
95!
96! 1241 2013-10-30 11:36:58Z heinze
97! output for profiles of ug and vg added
98! set surface_heatflux and surface_waterflux also in dependence on
99! large_scale_forcing
100! checks for nudging and large scale forcing from external file
101!
102! 1236 2013-09-27 07:21:13Z raasch
103! check number of spectra levels
104!
105! 1216 2013-08-26 09:31:42Z raasch
106! check for transpose_compute_overlap (temporary)
107!
108! 1214 2013-08-21 12:29:17Z kanani
109! additional check for simultaneous use of vertical grid stretching
110! and particle advection
111!
112! 1212 2013-08-15 08:46:27Z raasch
113! checks for poisfft_hybrid removed
114!
115! 1210 2013-08-14 10:58:20Z raasch
116! check for fftw
117!
118! 1179 2013-06-14 05:57:58Z raasch
119! checks and settings of buoyancy parameters and switches revised,
120! initial profile for rho added to hom (id=77)
121!
122! 1174 2013-05-31 10:28:08Z gryschka
123! Bugfix in computing initial profiles for ug, vg, lad, q in case of Atmosphere
124!
125! 1159 2013-05-21 11:58:22Z fricke
126! bc_lr/ns_dirneu/neudir removed
127!
128! 1115 2013-03-26 18:16:16Z hoffmann
129! unused variables removed
130! drizzle can be used without precipitation
131!
132! 1111 2013-03-08 23:54:10Z raasch
133! ibc_p_b = 2 removed
134!
135! 1103 2013-02-20 02:15:53Z raasch
136! Bugfix: turbulent inflow must not require cyclic fill in restart runs
137!
138! 1092 2013-02-02 11:24:22Z raasch
139! unused variables removed
140!
141! 1069 2012-11-28 16:18:43Z maronga
142! allow usage of topography in combination with cloud physics
143!
144! 1065 2012-11-22 17:42:36Z hoffmann
145! Bugfix: It is not allowed to use cloud_scheme = seifert_beheng without
146!         precipitation in order to save computational resources.
147!
148! 1060 2012-11-21 07:19:51Z raasch
149! additional check for parameter turbulent_inflow
150!
151! 1053 2012-11-13 17:11:03Z hoffmann
152! necessary changes for the new two-moment cloud physics scheme added:
153! - check cloud physics scheme (Kessler or Seifert and Beheng)
154! - plant_canopy is not allowed
155! - currently, only cache loop_optimization is allowed
156! - initial profiles of nr, qr
157! - boundary condition of nr, qr
158! - check output quantities (qr, nr, prr)
159!
160! 1036 2012-10-22 13:43:42Z raasch
161! code put under GPL (PALM 3.9)
162!
163! 1031/1034 2012-10-22 11:32:49Z raasch
164! check of netcdf4 parallel file support
165!
166! 1019 2012-09-28 06:46:45Z raasch
167! non-optimized version of prognostic_equations not allowed any more
168!
169! 1015 2012-09-27 09:23:24Z raasch
170! acc allowed for loop optimization,
171! checks for adjustment of mixing length to the Prandtl mixing length removed
172!
173! 1003 2012-09-14 14:35:53Z raasch
174! checks for cases with unequal subdomain sizes removed
175!
176! 1001 2012-09-13 14:08:46Z raasch
177! all actions concerning leapfrog- and upstream-spline-scheme removed
178!
179! 996 2012-09-07 10:41:47Z raasch
180! little reformatting
181
182! 978 2012-08-09 08:28:32Z fricke
183! setting of bc_lr/ns_dirneu/neudir
184! outflow damping layer removed
185! check for z0h*
186! check for pt_damping_width
187!
188! 964 2012-07-26 09:14:24Z raasch
189! check of old profil-parameters removed
190!
191! 940 2012-07-09 14:31:00Z raasch
192! checks for parameter neutral
193!
194! 924 2012-06-06 07:44:41Z maronga
195! Bugfix: preprocessor directives caused error during compilation
196!
197! 892 2012-05-02 13:51:44Z maronga
198! Bugfix for parameter file check ( excluding __netcdf4 )
199!
200! 866 2012-03-28 06:44:41Z raasch
201! use only 60% of the geostrophic wind as translation speed in case of Galilean
202! transformation and use_ug_for_galilei_tr = .T. in order to mimimize the
203! timestep
204!
205! 861 2012-03-26 14:18:34Z suehring
206! Check for topography and ws-scheme removed.
207! Check for loop_optimization = 'vector' and ws-scheme removed.
208!
209! 845 2012-03-07 10:23:05Z maronga
210! Bugfix: exclude __netcdf4 directive part from namelist file check compilation
211!
212! 828 2012-02-21 12:00:36Z raasch
213! check of collision_kernel extended
214!
215! 825 2012-02-19 03:03:44Z raasch
216! check for collision_kernel and curvature_solution_effects
217!
218! 809 2012-01-30 13:32:58Z maronga
219! Bugfix: replaced .AND. and .NOT. with && and ! in the preprocessor directives
220!
221! 807 2012-01-25 11:53:51Z maronga
222! New cpp directive "__check" implemented which is used by check_namelist_files
223!
224! Revision 1.1  1997/08/26 06:29:23  raasch
225! Initial revision
226!
227!
228! Description:
229! ------------
230! Check control parameters and deduce further quantities.
231!------------------------------------------------------------------------------!
232
233    USE arrays_3d
234    USE cloud_parameters
235    USE constants
236    USE control_parameters
237    USE dvrp_variables
238    USE grid_variables
239    USE indices
240    USE kinds
241    USE model_1d
242    USE netcdf_control
243    USE particle_attributes
244    USE pegrid
245    USE profil_parameter
246    USE spectrum
247    USE statistics
248    USE subsidence_mod
249    USE statistics
250    USE transpose_indices
251
252    IMPLICIT NONE
253
254    CHARACTER (LEN=1)   ::  sq                       !:
255    CHARACTER (LEN=6)   ::  var                      !:
256    CHARACTER (LEN=7)   ::  unit                     !:
257    CHARACTER (LEN=8)   ::  date                     !:
258    CHARACTER (LEN=10)  ::  time                     !:
259    CHARACTER (LEN=40)  ::  coupling_string          !:
260    CHARACTER (LEN=100) ::  action                   !:
261
262    INTEGER(iwp) ::  i                               !:
263    INTEGER(iwp) ::  ilen                            !:
264    INTEGER(iwp) ::  iremote = 0                     !:
265    INTEGER(iwp) ::  j                               !:
266    INTEGER(iwp) ::  k                               !:
267    INTEGER(iwp) ::  kk                              !:
268    INTEGER(iwp) ::  netcdf_data_format_save         !:
269    INTEGER(iwp) ::  position                        !:
270    INTEGER(iwp) ::  prec                            !:
271   
272    LOGICAL     ::  found                            !:
273    LOGICAL     ::  ldum                             !:
274   
275    REAL(wp)    ::  gradient                         !:
276    REAL(wp)    ::  remote = 0.0_wp                  !:
277    REAL(wp)    ::  simulation_time_since_reference  !:
278
279
280    CALL location_message( 'checking parameters', .FALSE. )
281
282!
283!-- Check for overlap combinations, which are not realized yet
284    IF ( transpose_compute_overlap )  THEN
285       IF ( numprocs == 1 )  STOP '+++ transpose-compute-overlap not implemented for single PE runs'
286#if defined( __openacc )
287       STOP '+++ transpose-compute-overlap not implemented for GPU usage'
288#endif
289    ENDIF
290
291!
292!-- Warning, if host is not set
293    IF ( host(1:1) == ' ' )  THEN
294       message_string = '"host" is not set. Please check that environment ' // &
295                        'variable "localhost" & is set before running PALM'
296       CALL message( 'check_parameters', 'PA0001', 0, 0, 0, 6, 0 )
297    ENDIF
298
299!
300!-- Check the coupling mode
301    IF ( coupling_mode /= 'uncoupled'            .AND.  &
302         coupling_mode /= 'atmosphere_to_ocean'  .AND.  &
303         coupling_mode /= 'ocean_to_atmosphere' )  THEN
304       message_string = 'illegal coupling mode: ' // TRIM( coupling_mode )
305       CALL message( 'check_parameters', 'PA0002', 1, 2, 0, 6, 0 )
306    ENDIF
307
308!
309!-- Check dt_coupling, restart_time, dt_restart, end_time, dx, dy, nx and ny
310    IF ( coupling_mode /= 'uncoupled')  THEN
311
312       IF ( dt_coupling == 9999999.9_wp )  THEN
313          message_string = 'dt_coupling is not set but required for coup' // &
314                           'ling mode "' //  TRIM( coupling_mode ) // '"'
315          CALL message( 'check_parameters', 'PA0003', 1, 2, 0, 6, 0 )
316       ENDIF
317
318#if defined( __parallel )
319
320!
321!--    NOTE: coupled runs have not been implemented in the check_namelist_files
322!--    program.
323!--    check_namelist_files will need the following information of the other
324!--    model (atmosphere/ocean).
325!       dt_coupling = remote
326!       dt_max = remote
327!       restart_time = remote
328!       dt_restart= remote
329!       simulation_time_since_reference = remote
330!       dx = remote
331
332
333#if ! defined( __check )
334       IF ( myid == 0 ) THEN
335          CALL MPI_SEND( dt_coupling, 1, MPI_REAL, target_id, 11, comm_inter, &
336                         ierr )
337          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 11, comm_inter, &
338                         status, ierr )
339       ENDIF
340       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
341#endif     
342       IF ( dt_coupling /= remote )  THEN
343          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
344                 '": dt_coupling = ', dt_coupling, '& is not equal to ',       &
345                 'dt_coupling_remote = ', remote
346          CALL message( 'check_parameters', 'PA0004', 1, 2, 0, 6, 0 )
347       ENDIF
348       IF ( dt_coupling <= 0.0_wp )  THEN
349#if ! defined( __check )
350          IF ( myid == 0  ) THEN
351             CALL MPI_SEND( dt_max, 1, MPI_REAL, target_id, 19, comm_inter, ierr )
352             CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 19, comm_inter, &
353                            status, ierr )
354          ENDIF   
355          CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
356#endif         
357          dt_coupling = MAX( dt_max, remote )
358          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
359                 '": dt_coupling <= 0.0 & is not allowed and is reset to ',    &
360                 'MAX(dt_max(A,O)) = ', dt_coupling
361          CALL message( 'check_parameters', 'PA0005', 0, 1, 0, 6, 0 )
362       ENDIF
363#if ! defined( __check )
364       IF ( myid == 0 ) THEN
365          CALL MPI_SEND( restart_time, 1, MPI_REAL, target_id, 12, comm_inter, &
366                         ierr )
367          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 12, comm_inter, &
368                         status, ierr )
369       ENDIF
370       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
371#endif     
372       IF ( restart_time /= remote )  THEN
373          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
374                 '": restart_time = ', restart_time, '& is not equal to ',     &
375                 'restart_time_remote = ', remote
376          CALL message( 'check_parameters', 'PA0006', 1, 2, 0, 6, 0 )
377       ENDIF
378#if ! defined( __check )
379       IF ( myid == 0 ) THEN
380          CALL MPI_SEND( dt_restart, 1, MPI_REAL, target_id, 13, comm_inter, &
381                         ierr )
382          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 13, comm_inter, &
383                         status, ierr )
384       ENDIF   
385       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
386#endif     
387       IF ( dt_restart /= remote )  THEN
388          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
389                 '": dt_restart = ', dt_restart, '& is not equal to ',         &
390                 'dt_restart_remote = ', remote
391          CALL message( 'check_parameters', 'PA0007', 1, 2, 0, 6, 0 )
392       ENDIF
393
394       simulation_time_since_reference = end_time - coupling_start_time
395#if ! defined( __check )
396       IF  ( myid == 0 ) THEN
397          CALL MPI_SEND( simulation_time_since_reference, 1, MPI_REAL, target_id, &
398                         14, comm_inter, ierr )
399          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 14, comm_inter, &
400                         status, ierr )   
401       ENDIF
402       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
403#endif     
404       IF ( simulation_time_since_reference /= remote )  THEN
405          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
406                 '": simulation_time_since_reference = ',                      &
407                 simulation_time_since_reference, '& is not equal to ',        &
408                 'simulation_time_since_reference_remote = ', remote
409          CALL message( 'check_parameters', 'PA0008', 1, 2, 0, 6, 0 )
410       ENDIF
411
412#if ! defined( __check )
413       IF ( myid == 0 ) THEN
414          CALL MPI_SEND( dx, 1, MPI_REAL, target_id, 15, comm_inter, ierr )
415          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 15, comm_inter, &
416                                                             status, ierr )
417       ENDIF
418       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
419
420#endif
421       IF ( coupling_mode == 'atmosphere_to_ocean') THEN
422
423          IF ( dx < remote ) THEN
424             WRITE( message_string, * ) 'coupling mode "', &
425                   TRIM( coupling_mode ),                  &
426           '": dx in Atmosphere is not equal to or not larger then dx in ocean'
427             CALL message( 'check_parameters', 'PA0009', 1, 2, 0, 6, 0 )
428          ENDIF
429
430          IF ( (nx_a+1)*dx /= (nx_o+1)*remote )  THEN
431             WRITE( message_string, * ) 'coupling mode "', &
432                    TRIM( coupling_mode ), &
433             '": Domain size in x-direction is not equal in ocean and atmosphere'
434             CALL message( 'check_parameters', 'PA0010', 1, 2, 0, 6, 0 )
435          ENDIF
436
437       ENDIF
438
439#if ! defined( __check )
440       IF ( myid == 0) THEN
441          CALL MPI_SEND( dy, 1, MPI_REAL, target_id, 16, comm_inter, ierr )
442          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 16, comm_inter, &
443                         status, ierr )
444       ENDIF
445       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
446#endif
447       IF ( coupling_mode == 'atmosphere_to_ocean') THEN
448
449          IF ( dy < remote )  THEN
450             WRITE( message_string, * ) 'coupling mode "', &
451                    TRIM( coupling_mode ), &
452                 '": dy in Atmosphere is not equal to or not larger then dy in ocean'
453             CALL message( 'check_parameters', 'PA0011', 1, 2, 0, 6, 0 )
454          ENDIF
455
456          IF ( (ny_a+1)*dy /= (ny_o+1)*remote )  THEN
457             WRITE( message_string, * ) 'coupling mode "', &
458                   TRIM( coupling_mode ), &
459             '": Domain size in y-direction is not equal in ocean and atmosphere'
460             CALL message( 'check_parameters', 'PA0012', 1, 2, 0, 6, 0 )
461          ENDIF
462
463          IF ( MOD(nx_o+1,nx_a+1) /= 0 )  THEN
464             WRITE( message_string, * ) 'coupling mode "', &
465                   TRIM( coupling_mode ), &
466             '": nx+1 in ocean is not divisible without remainder with nx+1 in', & 
467             ' atmosphere'
468             CALL message( 'check_parameters', 'PA0339', 1, 2, 0, 6, 0 )
469          ENDIF
470
471          IF ( MOD(ny_o+1,ny_a+1) /= 0 )  THEN
472             WRITE( message_string, * ) 'coupling mode "', &
473                   TRIM( coupling_mode ), &
474             '": ny+1 in ocean is not divisible without remainder with ny+1 in', & 
475             ' atmosphere'
476             CALL message( 'check_parameters', 'PA0340', 1, 2, 0, 6, 0 )
477          ENDIF
478
479       ENDIF
480#else
481       WRITE( message_string, * ) 'coupling requires PALM to be called with', &
482            ' ''mrun -K parallel'''
483       CALL message( 'check_parameters', 'PA0141', 1, 2, 0, 6, 0 )
484#endif
485    ENDIF
486
487#if defined( __parallel ) && ! defined ( __check )
488!
489!-- Exchange via intercommunicator
490    IF ( coupling_mode == 'atmosphere_to_ocean' .AND. myid == 0 )  THEN
491       CALL MPI_SEND( humidity, 1, MPI_LOGICAL, target_id, 19, comm_inter, &
492                      ierr )
493    ELSEIF ( coupling_mode == 'ocean_to_atmosphere' .AND. myid == 0)  THEN
494       CALL MPI_RECV( humidity_remote, 1, MPI_LOGICAL, target_id, 19, &
495                      comm_inter, status, ierr )
496    ENDIF
497    CALL MPI_BCAST( humidity_remote, 1, MPI_LOGICAL, 0, comm2d, ierr)
498   
499#endif
500
501
502!
503!-- Generate the file header which is used as a header for most of PALM's
504!-- output files
505    CALL DATE_AND_TIME( date, time )
506    run_date = date(7:8)//'-'//date(5:6)//'-'//date(3:4)
507    run_time = time(1:2)//':'//time(3:4)//':'//time(5:6)
508    IF ( coupling_mode == 'uncoupled' )  THEN
509       coupling_string = ''
510    ELSEIF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
511       coupling_string = ' coupled (atmosphere)'
512    ELSEIF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
513       coupling_string = ' coupled (ocean)'
514    ENDIF       
515
516    IF ( ensemble_member_nr /= 0 )  THEN
517       WRITE ( run_description_header,                                         &
518                  '(A,2X,A,2X,A,A,A,I2.2,A,2X,A,I2.2,2X,A,A,2X,A,1X,A)' )      &
519              TRIM( version ), TRIM( revision ), 'run: ',                      &
520              TRIM( run_identifier ), '.', runnr, TRIM( coupling_string ),     &
521              'en-no: ', ensemble_member_nr,'host: ', TRIM( host ),            &
522              run_date, run_time
523    ELSE
524       WRITE ( run_description_header,                                         &
525                  '(A,2X,A,2X,A,A,A,I2.2,A,2X,A,A,2X,A,1X,A)' )                &
526              TRIM( version ), TRIM( revision ), 'run: ',                      &
527              TRIM( run_identifier ), '.', runnr, TRIM( coupling_string ),     &
528              'host: ', TRIM( host ), run_date, run_time
529    ENDIF
530!
531!-- Check the general loop optimization method
532    IF ( loop_optimization == 'default' )  THEN
533       IF ( host(1:3) == 'nec' )  THEN
534          loop_optimization = 'vector'
535       ELSE
536          loop_optimization = 'cache'
537       ENDIF
538    ENDIF
539
540    SELECT CASE ( TRIM( loop_optimization ) )
541
542       CASE ( 'acc', 'cache', 'vector' )
543          CONTINUE
544
545       CASE DEFAULT
546          message_string = 'illegal value given for loop_optimization: "' // &
547                           TRIM( loop_optimization ) // '"'
548          CALL message( 'check_parameters', 'PA0013', 1, 2, 0, 6, 0 )
549
550    END SELECT
551
552!
553!-- Check if vertical grid stretching is used together with particles
554    IF ( dz_stretch_level < 100000.0_wp .AND. particle_advection )  THEN
555       message_string = 'Vertical grid stretching is not allowed together ' // &
556                        'with particle advection.'
557       CALL message( 'check_parameters', 'PA0017', 1, 2, 0, 6, 0 )
558    ENDIF
559
560!
561!--
562    IF ( use_particle_tails )  THEN
563       message_string = 'Particle tails are currently not available due ' //   &
564                        'to the new particle structure.'
565       CALL message( 'check_parameters', 'PA0392', 1, 2, 0, 6, 0 )
566    ENDIF
567
568!
569!-- Check topography setting (check for illegal parameter combinations)
570    IF ( topography /= 'flat' )  THEN
571       action = ' '
572       IF ( scalar_advec /= 'pw-scheme' .AND. scalar_advec /= 'ws-scheme')  THEN
573          WRITE( action, '(A,A)' )  'scalar_advec = ', scalar_advec
574       ENDIF
575       IF ( momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ws-scheme' ) &
576       THEN
577          WRITE( action, '(A,A)' )  'momentum_advec = ', momentum_advec
578       ENDIF
579       IF ( psolver == 'sor' )  THEN
580          WRITE( action, '(A,A)' )  'psolver = ', psolver
581       ENDIF
582       IF ( sloping_surface )  THEN
583          WRITE( action, '(A)' )  'sloping surface = .TRUE.'
584       ENDIF
585       IF ( galilei_transformation )  THEN
586          WRITE( action, '(A)' )  'galilei_transformation = .TRUE.'
587       ENDIF
588       IF ( cloud_physics )  THEN
589          WRITE( action, '(A)' )  'cloud_physics = .TRUE.'
590       ENDIF
591       IF ( cloud_droplets )  THEN
592          WRITE( action, '(A)' )  'cloud_droplets = .TRUE.'
593       ENDIF
594       IF ( .NOT. prandtl_layer )  THEN
595          WRITE( action, '(A)' )  'prandtl_layer = .FALSE.'
596       ENDIF
597       IF ( action /= ' ' )  THEN
598          message_string = 'a non-flat topography does not allow ' // &
599                           TRIM( action )
600          CALL message( 'check_parameters', 'PA0014', 1, 2, 0, 6, 0 )
601       ENDIF
602!
603!--    In case of non-flat topography, check whether the convention how to
604!--    define the topography grid has been set correctly, or whether the default
605!--    is applicable. If this is not possible, abort.
606       IF ( TRIM( topography_grid_convention ) == ' ' )  THEN
607          IF ( TRIM( topography ) /= 'single_building' .AND.  &
608               TRIM( topography ) /= 'single_street_canyon' .AND.  &
609               TRIM( topography ) /= 'read_from_file' )  THEN
610!--          The default value is not applicable here, because it is only valid
611!--          for the two standard cases 'single_building' and 'read_from_file'
612!--          defined in init_grid.
613             WRITE( message_string, * )  &
614                  'The value for "topography_grid_convention" ',  &
615                  'is not set. Its default value is & only valid for ',  &
616                  '"topography" = ''single_building'', ',  &
617                  '''single_street_canyon'' & or ''read_from_file''.',  &
618                  ' & Choose ''cell_edge'' or ''cell_center''.'
619             CALL message( 'user_check_parameters', 'PA0239', 1, 2, 0, 6, 0 )
620          ELSE
621!--          The default value is applicable here.
622!--          Set convention according to topography.
623             IF ( TRIM( topography ) == 'single_building' .OR.  &
624                  TRIM( topography ) == 'single_street_canyon' )  THEN
625                topography_grid_convention = 'cell_edge'
626             ELSEIF ( TRIM( topography ) == 'read_from_file' )  THEN
627                topography_grid_convention = 'cell_center'
628             ENDIF
629          ENDIF
630       ELSEIF ( TRIM( topography_grid_convention ) /= 'cell_edge' .AND.  &
631                TRIM( topography_grid_convention ) /= 'cell_center' )  THEN
632          WRITE( message_string, * )  &
633               'The value for "topography_grid_convention" is ', &
634               'not recognized. & Choose ''cell_edge'' or ''cell_center''.'
635          CALL message( 'user_check_parameters', 'PA0240', 1, 2, 0, 6, 0 )
636       ENDIF
637
638    ENDIF
639
640!
641!-- Check ocean setting
642    IF ( ocean )  THEN
643
644       action = ' '
645       IF ( action /= ' ' )  THEN
646          message_string = 'ocean = .T. does not allow ' // TRIM( action )
647          CALL message( 'check_parameters', 'PA0015', 1, 2, 0, 6, 0 )
648       ENDIF
649
650    ELSEIF ( TRIM( coupling_mode ) == 'uncoupled'  .AND.  &
651             TRIM( coupling_char ) == '_O' )  THEN
652
653!
654!--    Check whether an (uncoupled) atmospheric run has been declared as an
655!--    ocean run (this setting is done via mrun-option -y)
656
657       message_string = 'ocean = .F. does not allow coupling_char = "' // &
658                        TRIM( coupling_char ) // '" set by mrun-option "-y"'
659       CALL message( 'check_parameters', 'PA0317', 1, 2, 0, 6, 0 )
660
661    ENDIF
662!
663!-- Check cloud scheme
664    IF ( cloud_scheme == 'seifert_beheng' )  THEN
665       icloud_scheme = 0
666    ELSEIF ( cloud_scheme == 'kessler' )  THEN
667       icloud_scheme = 1
668    ELSE
669       message_string = 'unknown cloud microphysics scheme cloud_scheme ="' // &
670                        TRIM( cloud_scheme ) // '"'
671       CALL message( 'check_parameters', 'PA0357', 1, 2, 0, 6, 0 )
672    ENDIF
673!
674!-- Check whether there are any illegal values
675!-- Pressure solver:
676    IF ( psolver /= 'poisfft'  .AND. &
677         psolver /= 'sor'  .AND.  psolver /= 'multigrid' )  THEN
678       message_string = 'unknown solver for perturbation pressure: psolver' // &
679                        ' = "' // TRIM( psolver ) // '"'
680       CALL message( 'check_parameters', 'PA0016', 1, 2, 0, 6, 0 )
681    ENDIF
682
683    IF ( psolver == 'multigrid' )  THEN
684       IF ( cycle_mg == 'w' )  THEN
685          gamma_mg = 2
686       ELSEIF ( cycle_mg == 'v' )  THEN
687          gamma_mg = 1
688       ELSE
689          message_string = 'unknown multigrid cycle: cycle_mg = "' // &
690                           TRIM( cycle_mg ) // '"'
691          CALL message( 'check_parameters', 'PA0020', 1, 2, 0, 6, 0 )
692       ENDIF
693    ENDIF
694
695    IF ( fft_method /= 'singleton-algorithm'  .AND.  &
696         fft_method /= 'temperton-algorithm'  .AND.  &
697         fft_method /= 'fftw'                 .AND.  &
698         fft_method /= 'system-specific' )  THEN
699       message_string = 'unknown fft-algorithm: fft_method = "' // &
700                        TRIM( fft_method ) // '"'
701       CALL message( 'check_parameters', 'PA0021', 1, 2, 0, 6, 0 )
702    ENDIF
703   
704    IF( momentum_advec == 'ws-scheme' .AND. & 
705        .NOT. call_psolver_at_all_substeps  ) THEN
706        message_string = 'psolver must be called at each RK3 substep when "'//&
707                      TRIM(momentum_advec) // ' "is used for momentum_advec'
708        CALL message( 'check_parameters', 'PA0344', 1, 2, 0, 6, 0 )
709    END IF
710!
711!-- Advection schemes:
712    IF ( momentum_advec /= 'pw-scheme'  .AND.  momentum_advec /= 'ws-scheme' ) &
713    THEN
714       message_string = 'unknown advection scheme: momentum_advec = "' // &
715                        TRIM( momentum_advec ) // '"'
716       CALL message( 'check_parameters', 'PA0022', 1, 2, 0, 6, 0 )
717    ENDIF
718    IF ( ( momentum_advec == 'ws-scheme' .OR.  scalar_advec == 'ws-scheme' )   &
719           .AND. ( timestep_scheme == 'euler' .OR.                             &
720                   timestep_scheme == 'runge-kutta-2' ) )                      &
721    THEN
722       message_string = 'momentum_advec or scalar_advec = "' &
723         // TRIM( momentum_advec ) // '" is not allowed with timestep_scheme = "' // &
724         TRIM( timestep_scheme ) // '"'
725       CALL message( 'check_parameters', 'PA0023', 1, 2, 0, 6, 0 )
726    ENDIF
727    IF ( scalar_advec /= 'pw-scheme'  .AND.  scalar_advec /= 'ws-scheme' .AND. &
728         scalar_advec /= 'bc-scheme' )                                         &
729    THEN
730       message_string = 'unknown advection scheme: scalar_advec = "' // &
731                        TRIM( scalar_advec ) // '"'
732       CALL message( 'check_parameters', 'PA0024', 1, 2, 0, 6, 0 )
733    ENDIF
734    IF ( scalar_advec == 'bc-scheme'  .AND.  loop_optimization == 'cache' )    &
735    THEN
736       message_string = 'advection_scheme scalar_advec = "' &
737         // TRIM( scalar_advec ) // '" not implemented for & loop_optimization = "' // &
738         TRIM( loop_optimization ) // '"'
739       CALL message( 'check_parameters', 'PA0026', 1, 2, 0, 6, 0 )
740    ENDIF
741
742    IF ( use_sgs_for_particles  .AND.  .NOT. use_upstream_for_tke .AND.        &
743         scalar_advec /= 'ws-scheme' )  THEN
744       use_upstream_for_tke = .TRUE.
745       message_string = 'use_upstream_for_tke set .TRUE. because ' //          &
746                        'use_sgs_for_particles = .TRUE. '          //          &
747                        'and scalar_advec /= ws-scheme'
748       CALL message( 'check_parameters', 'PA0025', 0, 1, 0, 6, 0 )
749    ENDIF
750
751    IF ( use_sgs_for_particles  .AND.  curvature_solution_effects )  THEN
752       message_string = 'use_sgs_for_particles = .TRUE. not allowed with ' //  &
753                        'curvature_solution_effects = .TRUE.'
754       CALL message( 'check_parameters', 'PA0349', 1, 2, 0, 6, 0 )
755    ENDIF
756
757!
758!-- Set LOGICAL switches to enhance performance
759    IF ( momentum_advec == 'ws-scheme' )    ws_scheme_mom = .TRUE.
760    IF ( scalar_advec   == 'ws-scheme'   )  ws_scheme_sca = .TRUE.
761
762!
763!-- Timestep schemes:
764    SELECT CASE ( TRIM( timestep_scheme ) )
765
766       CASE ( 'euler' )
767          intermediate_timestep_count_max = 1
768
769       CASE ( 'runge-kutta-2' )
770          intermediate_timestep_count_max = 2
771
772       CASE ( 'runge-kutta-3' )
773          intermediate_timestep_count_max = 3
774
775       CASE DEFAULT
776          message_string = 'unknown timestep scheme: timestep_scheme = "' //   &
777                           TRIM( timestep_scheme ) // '"'
778          CALL message( 'check_parameters', 'PA0027', 1, 2, 0, 6, 0 )
779
780    END SELECT
781
782    IF ( (momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ws-scheme')   &
783         .AND. timestep_scheme(1:5) == 'runge' ) THEN
784       message_string = 'momentum advection scheme "' // &
785                        TRIM( momentum_advec ) // '" & does not work with ' // &
786                        'timestep_scheme "' // TRIM( timestep_scheme ) // '"'
787       CALL message( 'check_parameters', 'PA0029', 1, 2, 0, 6, 0 )
788    ENDIF
789
790!
791!-- Collision kernels:
792    SELECT CASE ( TRIM( collision_kernel ) )
793
794       CASE ( 'hall', 'hall_fast' )
795          hall_kernel = .TRUE.
796
797       CASE ( 'palm' )
798          palm_kernel = .TRUE.
799
800       CASE ( 'wang', 'wang_fast' )
801          wang_kernel = .TRUE.
802
803       CASE ( 'none' )
804
805
806       CASE DEFAULT
807          message_string = 'unknown collision kernel: collision_kernel = "' // &
808                           TRIM( collision_kernel ) // '"'
809          CALL message( 'check_parameters', 'PA0350', 1, 2, 0, 6, 0 )
810
811    END SELECT
812    IF ( collision_kernel(6:9) == 'fast' )  use_kernel_tables = .TRUE.
813
814    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.  &
815         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
816!
817!--    No restart run: several initialising actions are possible
818       action = initializing_actions
819       DO WHILE ( TRIM( action ) /= '' )
820          position = INDEX( action, ' ' )
821          SELECT CASE ( action(1:position-1) )
822
823             CASE ( 'set_constant_profiles', 'set_1d-model_profiles', &
824                    'by_user', 'initialize_vortex',     'initialize_ptanom' )
825                action = action(position+1:)
826
827             CASE DEFAULT
828                message_string = 'initializing_action = "' // &
829                                 TRIM( action ) // '" unkown or not allowed'
830                CALL message( 'check_parameters', 'PA0030', 1, 2, 0, 6, 0 )
831
832          END SELECT
833       ENDDO
834    ENDIF
835
836    IF ( TRIM( initializing_actions ) == 'initialize_vortex' .AND. &
837         conserve_volume_flow ) THEN
838         message_string = 'initializing_actions = "initialize_vortex"' // &
839                        ' ist not allowed with conserve_volume_flow = .T.'
840       CALL message( 'check_parameters', 'PA0343', 1, 2, 0, 6, 0 )
841    ENDIF       
842
843
844    IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0  .AND. &
845         INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
846       message_string = 'initializing_actions = "set_constant_profiles"' // &
847                        ' and "set_1d-model_profiles" are not allowed ' //  &
848                        'simultaneously'
849       CALL message( 'check_parameters', 'PA0031', 1, 2, 0, 6, 0 )
850    ENDIF
851
852    IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0  .AND. &
853         INDEX( initializing_actions, 'by_user' ) /= 0 )  THEN
854       message_string = 'initializing_actions = "set_constant_profiles"' // &
855                        ' and "by_user" are not allowed simultaneously'
856       CALL message( 'check_parameters', 'PA0032', 1, 2, 0, 6, 0 )
857    ENDIF
858
859    IF ( INDEX( initializing_actions, 'by_user' ) /= 0  .AND. &
860         INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
861       message_string = 'initializing_actions = "by_user" and ' // &
862                        '"set_1d-model_profiles" are not allowed simultaneously'
863       CALL message( 'check_parameters', 'PA0033', 1, 2, 0, 6, 0 )
864    ENDIF
865
866    IF ( cloud_physics  .AND.  .NOT. humidity )  THEN
867       WRITE( message_string, * ) 'cloud_physics = ', cloud_physics, ' is ', &
868              'not allowed with humidity = ', humidity
869       CALL message( 'check_parameters', 'PA0034', 1, 2, 0, 6, 0 )
870    ENDIF
871
872    IF ( precipitation  .AND.  .NOT.  cloud_physics )  THEN
873       WRITE( message_string, * ) 'precipitation = ', precipitation, ' is ', &
874              'not allowed with cloud_physics = ', cloud_physics
875       CALL message( 'check_parameters', 'PA0035', 1, 2, 0, 6, 0 )
876    ENDIF
877
878    IF ( humidity  .AND.  sloping_surface )  THEN
879       message_string = 'humidity = .TRUE. and sloping_surface = .TRUE. ' // &
880                        'are not allowed simultaneously'
881       CALL message( 'check_parameters', 'PA0036', 1, 2, 0, 6, 0 )
882    ENDIF
883
884    IF ( passive_scalar  .AND.  humidity )  THEN
885       message_string = 'humidity = .TRUE. and passive_scalar = .TRUE. ' // &
886                        'is not allowed simultaneously'
887       CALL message( 'check_parameters', 'PA0038', 1, 2, 0, 6, 0 )
888    ENDIF
889
890    IF ( plant_canopy .AND. ( drag_coefficient == 0.0_wp ) ) THEN
891       message_string = 'plant_canopy = .TRUE. requires a non-zero drag ' // &
892                        'coefficient & given value is drag_coefficient = 0.0'
893       CALL message( 'check_parameters', 'PA0041', 1, 2, 0, 6, 0 )
894    ENDIF
895
896    IF ( plant_canopy  .AND.  cloud_physics  .AND.  icloud_scheme == 0 ) THEN
897       message_string = 'plant_canopy = .TRUE. requires cloud_scheme /=' //  &
898                        ' seifert_beheng'
899       CALL message( 'check_parameters', 'PA0360', 1, 2, 0, 6, 0 )
900    ENDIF
901
902    IF ( .NOT. ( loop_optimization == 'cache'  .OR.                            &
903                 loop_optimization == 'vector' )                               &
904         .AND.  cloud_physics  .AND.  icloud_scheme == 0 )  THEN
905       message_string = 'cloud_scheme = seifert_beheng requires ' // &
906                        'loop_optimization = "cache" or "vector"'
907       CALL message( 'check_parameters', 'PA0362', 1, 2, 0, 6, 0 )
908    ENDIF 
909
910!
911!-- In case of no model continuation run, check initialising parameters and
912!-- deduce further quantities
913    IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
914
915!
916!--    Initial profiles for 1D and 3D model, respectively (u,v further below)
917       pt_init = pt_surface
918       IF ( humidity )  THEN
919          q_init  = q_surface
920       ENDIF
921       IF ( ocean )           sa_init = sa_surface
922       IF ( passive_scalar )  q_init  = s_surface
923       IF ( plant_canopy )    lad = 0.0_wp
924
925!
926!--
927!--    If required, compute initial profile of the geostrophic wind
928!--    (component ug)
929       i = 1
930       gradient = 0.0_wp
931
932       IF ( .NOT. ocean )  THEN
933
934          ug_vertical_gradient_level_ind(1) = 0
935          ug(0) = ug_surface
936          DO  k = 1, nzt+1
937             IF ( i < 11 ) THEN
938                IF ( ug_vertical_gradient_level(i) < zu(k)  .AND. &
939                     ug_vertical_gradient_level(i) >= 0.0_wp )  THEN
940                   gradient = ug_vertical_gradient(i) / 100.0_wp
941                   ug_vertical_gradient_level_ind(i) = k - 1
942                   i = i + 1
943                ENDIF
944             ENDIF       
945             IF ( gradient /= 0.0_wp )  THEN
946                IF ( k /= 1 )  THEN
947                   ug(k) = ug(k-1) + dzu(k) * gradient
948                ELSE
949                   ug(k) = ug_surface + dzu(k) * gradient
950                ENDIF
951             ELSE
952                ug(k) = ug(k-1)
953             ENDIF
954          ENDDO
955
956       ELSE
957
958          ug_vertical_gradient_level_ind(1) = nzt+1
959          ug(nzt+1) = ug_surface
960          DO  k = nzt, nzb, -1
961             IF ( i < 11 ) THEN
962                IF ( ug_vertical_gradient_level(i) > zu(k)  .AND. &
963                     ug_vertical_gradient_level(i) <= 0.0_wp )  THEN
964                   gradient = ug_vertical_gradient(i) / 100.0_wp
965                   ug_vertical_gradient_level_ind(i) = k + 1
966                   i = i + 1
967                ENDIF
968             ENDIF
969             IF ( gradient /= 0.0_wp )  THEN
970                IF ( k /= nzt )  THEN
971                   ug(k) = ug(k+1) - dzu(k+1) * gradient
972                ELSE
973                   ug(k)   = ug_surface - 0.5_wp * dzu(k+1) * gradient
974                   ug(k+1) = ug_surface + 0.5_wp * dzu(k+1) * gradient
975                ENDIF
976             ELSE
977                ug(k) = ug(k+1)
978             ENDIF
979          ENDDO
980
981       ENDIF
982
983!
984!--    In case of no given gradients for ug, choose a zero gradient
985       IF ( ug_vertical_gradient_level(1) == -9999999.9_wp )  THEN
986          ug_vertical_gradient_level(1) = 0.0_wp
987       ENDIF 
988
989!
990!--
991!--    If required, compute initial profile of the geostrophic wind
992!--    (component vg)
993       i = 1
994       gradient = 0.0_wp
995
996       IF ( .NOT. ocean )  THEN
997
998          vg_vertical_gradient_level_ind(1) = 0
999          vg(0) = vg_surface
1000          DO  k = 1, nzt+1
1001             IF ( i < 11 ) THEN
1002                IF ( vg_vertical_gradient_level(i) < zu(k)  .AND. &
1003                     vg_vertical_gradient_level(i) >= 0.0_wp )  THEN
1004                   gradient = vg_vertical_gradient(i) / 100.0_wp
1005                   vg_vertical_gradient_level_ind(i) = k - 1
1006                   i = i + 1
1007                ENDIF
1008             ENDIF
1009             IF ( gradient /= 0.0_wp )  THEN
1010                IF ( k /= 1 )  THEN
1011                   vg(k) = vg(k-1) + dzu(k) * gradient
1012                ELSE
1013                   vg(k) = vg_surface + dzu(k) * gradient
1014                ENDIF
1015             ELSE
1016                vg(k) = vg(k-1)
1017             ENDIF
1018          ENDDO
1019
1020       ELSE
1021
1022          vg_vertical_gradient_level_ind(1) = nzt+1
1023          vg(nzt+1) = vg_surface
1024          DO  k = nzt, nzb, -1
1025             IF ( i < 11 ) THEN
1026                IF ( vg_vertical_gradient_level(i) > zu(k)  .AND. &
1027                     vg_vertical_gradient_level(i) <= 0.0_wp )  THEN
1028                   gradient = vg_vertical_gradient(i) / 100.0_wp
1029                   vg_vertical_gradient_level_ind(i) = k + 1
1030                   i = i + 1
1031                ENDIF
1032             ENDIF
1033             IF ( gradient /= 0.0_wp )  THEN
1034                IF ( k /= nzt )  THEN
1035                   vg(k) = vg(k+1) - dzu(k+1) * gradient
1036                ELSE
1037                   vg(k)   = vg_surface - 0.5_wp * dzu(k+1) * gradient
1038                   vg(k+1) = vg_surface + 0.5_wp * dzu(k+1) * gradient
1039                ENDIF
1040             ELSE
1041                vg(k) = vg(k+1)
1042             ENDIF
1043          ENDDO
1044
1045       ENDIF
1046
1047!
1048!--    In case of no given gradients for vg, choose a zero gradient
1049       IF ( vg_vertical_gradient_level(1) == -9999999.9_wp )  THEN
1050          vg_vertical_gradient_level(1) = 0.0_wp
1051       ENDIF
1052
1053!
1054!--    Let the initial wind profiles be the calculated ug/vg profiles or
1055!--    interpolate them from wind profile data (if given)
1056       IF ( u_profile(1) == 9999999.9_wp  .AND.  v_profile(1) == 9999999.9_wp )  THEN
1057
1058          u_init = ug
1059          v_init = vg
1060
1061       ELSEIF ( u_profile(1) == 0.0_wp  .AND.  v_profile(1) == 0.0_wp )  THEN
1062
1063          IF ( uv_heights(1) /= 0.0_wp )  THEN
1064             message_string = 'uv_heights(1) must be 0.0'
1065             CALL message( 'check_parameters', 'PA0345', 1, 2, 0, 6, 0 )
1066          ENDIF
1067
1068          use_prescribed_profile_data = .TRUE.
1069
1070          kk = 1
1071          u_init(0) = 0.0_wp
1072          v_init(0) = 0.0_wp
1073
1074          DO  k = 1, nz+1
1075
1076             IF ( kk < 100 )  THEN
1077                DO WHILE ( uv_heights(kk+1) <= zu(k) )
1078                   kk = kk + 1
1079                   IF ( kk == 100 )  EXIT
1080                ENDDO
1081             ENDIF
1082
1083             IF ( kk < 100 .AND. uv_heights(kk+1) /= 9999999.9_wp )  THEN
1084                u_init(k) = u_profile(kk) + ( zu(k) - uv_heights(kk) ) /       &
1085                                       ( uv_heights(kk+1) - uv_heights(kk) ) * &
1086                                       ( u_profile(kk+1) - u_profile(kk) )
1087                v_init(k) = v_profile(kk) + ( zu(k) - uv_heights(kk) ) /       &
1088                                       ( uv_heights(kk+1) - uv_heights(kk) ) * &
1089                                       ( v_profile(kk+1) - v_profile(kk) )
1090             ELSE
1091                u_init(k) = u_profile(kk)
1092                v_init(k) = v_profile(kk)
1093             ENDIF
1094
1095          ENDDO
1096
1097       ELSE
1098
1099          message_string = 'u_profile(1) and v_profile(1) must be 0.0'
1100          CALL message( 'check_parameters', 'PA0346', 1, 2, 0, 6, 0 )
1101
1102       ENDIF
1103
1104!
1105!--    Compute initial temperature profile using the given temperature gradients
1106       IF ( .NOT. neutral )  THEN
1107
1108          i = 1
1109          gradient = 0.0_wp
1110
1111          IF ( .NOT. ocean )  THEN
1112
1113             pt_vertical_gradient_level_ind(1) = 0
1114             DO  k = 1, nzt+1
1115                IF ( i < 11 ) THEN
1116                   IF ( pt_vertical_gradient_level(i) < zu(k)  .AND. &
1117                        pt_vertical_gradient_level(i) >= 0.0_wp )  THEN
1118                      gradient = pt_vertical_gradient(i) / 100.0_wp
1119                      pt_vertical_gradient_level_ind(i) = k - 1
1120                      i = i + 1
1121                   ENDIF
1122                ENDIF
1123                IF ( gradient /= 0.0_wp )  THEN
1124                   IF ( k /= 1 )  THEN
1125                      pt_init(k) = pt_init(k-1) + dzu(k) * gradient
1126                   ELSE
1127                      pt_init(k) = pt_surface   + dzu(k) * gradient
1128                   ENDIF
1129                ELSE
1130                   pt_init(k) = pt_init(k-1)
1131                ENDIF
1132             ENDDO
1133
1134          ELSE
1135
1136             pt_vertical_gradient_level_ind(1) = nzt+1
1137             DO  k = nzt, 0, -1
1138                IF ( i < 11 ) THEN
1139                   IF ( pt_vertical_gradient_level(i) > zu(k)  .AND. &
1140                        pt_vertical_gradient_level(i) <= 0.0_wp )  THEN
1141                      gradient = pt_vertical_gradient(i) / 100.0_wp
1142                      pt_vertical_gradient_level_ind(i) = k + 1
1143                      i = i + 1
1144                   ENDIF
1145                ENDIF
1146                IF ( gradient /= 0.0_wp )  THEN
1147                   IF ( k /= nzt )  THEN
1148                      pt_init(k) = pt_init(k+1) - dzu(k+1) * gradient
1149                   ELSE
1150                      pt_init(k)   = pt_surface - 0.5_wp * dzu(k+1) * gradient
1151                      pt_init(k+1) = pt_surface + 0.5_wp * dzu(k+1) * gradient
1152                   ENDIF
1153                ELSE
1154                   pt_init(k) = pt_init(k+1)
1155                ENDIF
1156             ENDDO
1157
1158          ENDIF
1159
1160       ENDIF
1161
1162!
1163!--    In case of no given temperature gradients, choose gradient of neutral
1164!--    stratification
1165       IF ( pt_vertical_gradient_level(1) == -9999999.9_wp )  THEN
1166          pt_vertical_gradient_level(1) = 0.0_wp
1167       ENDIF
1168
1169!
1170!--    Store temperature gradient at the top boundary for possible Neumann
1171!--    boundary condition
1172       bc_pt_t_val = ( pt_init(nzt+1) - pt_init(nzt) ) / dzu(nzt+1)
1173
1174!
1175!--    If required, compute initial humidity or scalar profile using the given
1176!--    humidity/scalar gradient. In case of scalar transport, initially store
1177!--    values of the scalar parameters on humidity parameters
1178       IF ( passive_scalar )  THEN
1179          bc_q_b                    = bc_s_b
1180          bc_q_t                    = bc_s_t
1181          q_surface                 = s_surface
1182          q_surface_initial_change  = s_surface_initial_change
1183          q_vertical_gradient       = s_vertical_gradient
1184          q_vertical_gradient_level = s_vertical_gradient_level
1185          surface_waterflux         = surface_scalarflux
1186          wall_humidityflux         = wall_scalarflux
1187       ENDIF
1188
1189       IF ( humidity  .OR.  passive_scalar )  THEN
1190
1191          i = 1
1192          gradient = 0.0_wp
1193          q_vertical_gradient_level_ind(1) = 0
1194          DO  k = 1, nzt+1
1195             IF ( i < 11 ) THEN
1196                IF ( q_vertical_gradient_level(i) < zu(k)  .AND. &
1197                     q_vertical_gradient_level(i) >= 0.0_wp )  THEN
1198                   gradient = q_vertical_gradient(i) / 100.0_wp
1199                   q_vertical_gradient_level_ind(i) = k - 1
1200                   i = i + 1
1201                ENDIF
1202             ENDIF
1203             IF ( gradient /= 0.0_wp )  THEN
1204                IF ( k /= 1 )  THEN
1205                   q_init(k) = q_init(k-1) + dzu(k) * gradient
1206                ELSE
1207                   q_init(k) = q_init(k-1) + dzu(k) * gradient
1208                ENDIF
1209             ELSE
1210                q_init(k) = q_init(k-1)
1211             ENDIF
1212!
1213!--          Avoid negative humidities
1214             IF ( q_init(k) < 0.0_wp )  THEN
1215                q_init(k) = 0.0_wp
1216             ENDIF
1217          ENDDO
1218
1219!
1220!--       In case of no given humidity gradients, choose zero gradient
1221!--       conditions
1222          IF ( q_vertical_gradient_level(1) == -1.0_wp )  THEN
1223             q_vertical_gradient_level(1) = 0.0_wp
1224          ENDIF
1225!
1226!--       Store humidity, rain water content and rain drop concentration
1227!--       gradient at the top boundary for possile Neumann boundary condition
1228          bc_q_t_val  = ( q_init(nzt+1) - q_init(nzt) ) / dzu(nzt+1)
1229       ENDIF
1230
1231!
1232!--    If required, compute initial salinity profile using the given salinity
1233!--    gradients
1234       IF ( ocean )  THEN
1235
1236          i = 1
1237          gradient = 0.0_wp
1238
1239          sa_vertical_gradient_level_ind(1) = nzt+1
1240          DO  k = nzt, 0, -1
1241             IF ( i < 11 ) THEN
1242                IF ( sa_vertical_gradient_level(i) > zu(k)  .AND. &
1243                     sa_vertical_gradient_level(i) <= 0.0_wp )  THEN
1244                   gradient = sa_vertical_gradient(i) / 100.0_wp
1245                   sa_vertical_gradient_level_ind(i) = k + 1
1246                   i = i + 1
1247                ENDIF
1248             ENDIF
1249             IF ( gradient /= 0.0_wp )  THEN
1250                IF ( k /= nzt )  THEN
1251                   sa_init(k) = sa_init(k+1) - dzu(k+1) * gradient
1252                ELSE
1253                   sa_init(k)   = sa_surface - 0.5_wp * dzu(k+1) * gradient
1254                   sa_init(k+1) = sa_surface + 0.5_wp * dzu(k+1) * gradient
1255                ENDIF
1256             ELSE
1257                sa_init(k) = sa_init(k+1)
1258             ENDIF
1259          ENDDO
1260
1261       ENDIF
1262
1263!
1264!--    If required compute the profile of leaf area density used in the plant
1265!--    canopy model
1266       IF ( plant_canopy ) THEN
1267       
1268          i = 1
1269          gradient = 0.0_wp
1270
1271          IF ( .NOT. ocean ) THEN
1272
1273             lad(0) = lad_surface
1274 
1275             lad_vertical_gradient_level_ind(1) = 0
1276             DO k = 1, pch_index
1277                IF ( i < 11 ) THEN
1278                   IF ( lad_vertical_gradient_level(i) < zu(k) .AND.  &
1279                        lad_vertical_gradient_level(i) >= 0.0_wp ) THEN
1280                      gradient = lad_vertical_gradient(i)
1281                      lad_vertical_gradient_level_ind(i) = k - 1
1282                      i = i + 1
1283                   ENDIF
1284                ENDIF
1285                IF ( gradient /= 0.0_wp ) THEN
1286                   IF ( k /= 1 ) THEN
1287                      lad(k) = lad(k-1) + dzu(k) * gradient
1288                   ELSE
1289                      lad(k) = lad_surface + dzu(k) *gradient
1290                   ENDIF
1291                ELSE
1292                   lad(k) = lad(k-1)
1293                ENDIF
1294             ENDDO
1295
1296          ENDIF
1297
1298!
1299!--       In case of no given leaf area density gradients, choose a vanishing
1300!--       gradient
1301          IF ( lad_vertical_gradient_level(1) == -9999999.9_wp ) THEN
1302             lad_vertical_gradient_level(1) = 0.0_wp
1303          ENDIF
1304
1305       ENDIF
1306         
1307    ENDIF
1308
1309!
1310!-- Check if the control parameter use_subsidence_tendencies is used correctly
1311    IF ( use_subsidence_tendencies  .AND.  .NOT. large_scale_subsidence )  THEN
1312       message_string = 'The usage of use_subsidence_tendencies ' // &
1313                            'requires large_scale_subsidence = .T..'
1314       CALL message( 'check_parameters', 'PA0396', 1, 2, 0, 6, 0 )
1315    ELSEIF ( use_subsidence_tendencies  .AND.  .NOT. large_scale_forcing )  THEN
1316       message_string = 'The usage of use_subsidence_tendencies ' // &
1317                            'requires large_scale_forcing = .T..'
1318       CALL message( 'check_parameters', 'PA0397', 1, 2, 0, 6, 0 )
1319    ENDIF
1320
1321!
1322!-- Initialize large scale subsidence if required
1323    If ( large_scale_subsidence )  THEN
1324       IF ( subs_vertical_gradient_level(1) /= -9999999.9_wp .AND. &
1325                                     .NOT. large_scale_forcing )  THEN
1326          CALL init_w_subsidence
1327       ENDIF
1328!
1329!--    In case large_scale_forcing is used, profiles for subsidence velocity
1330!--    are read in from file LSF_DATA
1331
1332       IF ( subs_vertical_gradient_level(1) == -9999999.9_wp .AND. &
1333                                     .NOT. large_scale_forcing )  THEN
1334          message_string = 'There is no default large scale vertical ' // &
1335                           'velocity profile set. Specify the subsidence ' // &
1336                           'velocity profile via subs_vertical_gradient and ' // &
1337                           'subs_vertical_gradient_level.'
1338          CALL message( 'check_parameters', 'PA0380', 1, 2, 0, 6, 0 )
1339       ENDIF
1340    ELSE
1341        IF ( subs_vertical_gradient_level(1) /= -9999999.9_wp )  THEN
1342           message_string = 'Enable usage of large scale subsidence by ' // &
1343                            'setting large_scale_subsidence = .T..'
1344          CALL message( 'check_parameters', 'PA0381', 1, 2, 0, 6, 0 )
1345        ENDIF
1346    ENDIF   
1347
1348!
1349!-- Compute Coriolis parameter
1350    f  = 2.0_wp * omega * SIN( phi / 180.0_wp * pi )
1351    fs = 2.0_wp * omega * COS( phi / 180.0_wp * pi )
1352
1353!
1354!-- Check and set buoyancy related parameters and switches
1355    IF ( reference_state == 'horizontal_average' )  THEN
1356       CONTINUE
1357    ELSEIF ( reference_state == 'initial_profile' )  THEN
1358       use_initial_profile_as_reference = .TRUE.
1359    ELSEIF ( reference_state == 'single_value' )  THEN
1360       use_single_reference_value = .TRUE.
1361       IF ( pt_reference == 9999999.9_wp )  pt_reference = pt_surface
1362       vpt_reference = pt_reference * ( 1.0_wp + 0.61_wp * q_surface )
1363    ELSE
1364       message_string = 'illegal value for reference_state: "' // &
1365                        TRIM( reference_state ) // '"'
1366       CALL message( 'check_parameters', 'PA0056', 1, 2, 0, 6, 0 )
1367    ENDIF
1368
1369!
1370!-- Ocean runs always use reference values in the buoyancy term
1371    IF ( ocean )  THEN
1372       reference_state = 'single_value'
1373       use_single_reference_value = .TRUE.
1374    ENDIF
1375
1376!
1377!-- Sign of buoyancy/stability terms
1378    IF ( ocean )  atmos_ocean_sign = -1.0_wp
1379
1380!
1381!-- Ocean version must use flux boundary conditions at the top
1382    IF ( ocean .AND. .NOT. use_top_fluxes )  THEN
1383       message_string = 'use_top_fluxes must be .TRUE. in ocean mode'
1384       CALL message( 'check_parameters', 'PA0042', 1, 2, 0, 6, 0 )
1385    ENDIF
1386
1387!
1388!-- In case of a given slope, compute the relevant quantities
1389    IF ( alpha_surface /= 0.0_wp )  THEN
1390       IF ( ABS( alpha_surface ) > 90.0_wp )  THEN
1391          WRITE( message_string, * ) 'ABS( alpha_surface = ', alpha_surface, &
1392                                     ' ) must be < 90.0'
1393          CALL message( 'check_parameters', 'PA0043', 1, 2, 0, 6, 0 )
1394       ENDIF
1395       sloping_surface = .TRUE.
1396       cos_alpha_surface = COS( alpha_surface / 180.0_wp * pi )
1397       sin_alpha_surface = SIN( alpha_surface / 180.0_wp * pi )
1398    ENDIF
1399
1400!
1401!-- Check time step and cfl_factor
1402    IF ( dt /= -1.0_wp )  THEN
1403       IF ( dt <= 0.0_wp  .AND.  dt /= -1.0_wp )  THEN
1404          WRITE( message_string, * ) 'dt = ', dt , ' <= 0.0'
1405          CALL message( 'check_parameters', 'PA0044', 1, 2, 0, 6, 0 )
1406       ENDIF
1407       dt_3d = dt
1408       dt_fixed = .TRUE.
1409    ENDIF
1410
1411    IF ( cfl_factor <= 0.0_wp  .OR.  cfl_factor > 1.0_wp )  THEN
1412       IF ( cfl_factor == -1.0_wp )  THEN
1413          IF ( timestep_scheme == 'runge-kutta-2' )  THEN
1414             cfl_factor = 0.8_wp
1415          ELSEIF ( timestep_scheme == 'runge-kutta-3' )  THEN
1416             cfl_factor = 0.9_wp
1417          ELSE
1418             cfl_factor = 0.9_wp
1419          ENDIF
1420       ELSE
1421          WRITE( message_string, * ) 'cfl_factor = ', cfl_factor, &
1422                 ' out of range & 0.0 < cfl_factor <= 1.0 is required'
1423          CALL message( 'check_parameters', 'PA0045', 1, 2, 0, 6, 0 )
1424       ENDIF
1425    ENDIF
1426
1427!
1428!-- Store simulated time at begin
1429    simulated_time_at_begin = simulated_time
1430
1431!
1432!-- Store reference time for coupled runs and change the coupling flag,
1433!-- if ...
1434    IF ( simulated_time == 0.0_wp )  THEN
1435       IF ( coupling_start_time == 0.0_wp )  THEN
1436          time_since_reference_point = 0.0_wp
1437       ELSEIF ( time_since_reference_point < 0.0_wp )  THEN
1438          run_coupled = .FALSE.
1439       ENDIF
1440    ENDIF
1441
1442!
1443!-- Set wind speed in the Galilei-transformed system
1444    IF ( galilei_transformation )  THEN
1445       IF ( use_ug_for_galilei_tr .AND.                     &
1446            ug_vertical_gradient_level(1) == 0.0_wp  .AND.  &
1447            ug_vertical_gradient(1) == 0.0_wp  .AND.        & 
1448            vg_vertical_gradient_level(1) == 0.0_wp  .AND.  &
1449            vg_vertical_gradient(1) == 0.0_wp )  THEN
1450          u_gtrans = ug_surface * 0.6_wp
1451          v_gtrans = vg_surface * 0.6_wp
1452       ELSEIF ( use_ug_for_galilei_tr  .AND.                     &
1453                ( ug_vertical_gradient_level(1) /= 0.0_wp  .OR.  &
1454                ug_vertical_gradient(1) /= 0.0_wp ) )  THEN
1455          message_string = 'baroclinity (ug) not allowed simultaneously' // &
1456                           ' with galilei transformation'
1457          CALL message( 'check_parameters', 'PA0046', 1, 2, 0, 6, 0 )
1458       ELSEIF ( use_ug_for_galilei_tr  .AND.                     &
1459                ( vg_vertical_gradient_level(1) /= 0.0_wp  .OR.  &
1460                vg_vertical_gradient(1) /= 0.0_wp ) )  THEN
1461          message_string = 'baroclinity (vg) not allowed simultaneously' // &
1462                           ' with galilei transformation'
1463          CALL message( 'check_parameters', 'PA0047', 1, 2, 0, 6, 0 )
1464       ELSE
1465          message_string = 'variable translation speed used for galilei-' // &
1466             'transformation, which may cause & instabilities in stably ' // &
1467             'stratified regions'
1468          CALL message( 'check_parameters', 'PA0048', 0, 1, 0, 6, 0 )
1469       ENDIF
1470    ENDIF
1471
1472!
1473!-- In case of using a prandtl-layer, calculated (or prescribed) surface
1474!-- fluxes have to be used in the diffusion-terms
1475    IF ( prandtl_layer )  use_surface_fluxes = .TRUE.
1476
1477!
1478!-- Check boundary conditions and set internal variables:
1479!-- Lateral boundary conditions
1480    IF ( bc_lr /= 'cyclic'  .AND.  bc_lr /= 'dirichlet/radiation'  .AND. &
1481         bc_lr /= 'radiation/dirichlet' )  THEN
1482       message_string = 'unknown boundary condition: bc_lr = "' // &
1483                        TRIM( bc_lr ) // '"'
1484       CALL message( 'check_parameters', 'PA0049', 1, 2, 0, 6, 0 )
1485    ENDIF
1486    IF ( bc_ns /= 'cyclic'  .AND.  bc_ns /= 'dirichlet/radiation'  .AND. &
1487         bc_ns /= 'radiation/dirichlet' )  THEN
1488       message_string = 'unknown boundary condition: bc_ns = "' // &
1489                        TRIM( bc_ns ) // '"'
1490       CALL message( 'check_parameters', 'PA0050', 1, 2, 0, 6, 0 )
1491    ENDIF
1492
1493!
1494!-- Internal variables used for speed optimization in if clauses
1495    IF ( bc_lr /= 'cyclic' )               bc_lr_cyc    = .FALSE.
1496    IF ( bc_lr == 'dirichlet/radiation' )  bc_lr_dirrad = .TRUE.
1497    IF ( bc_lr == 'radiation/dirichlet' )  bc_lr_raddir = .TRUE.
1498    IF ( bc_ns /= 'cyclic' )               bc_ns_cyc    = .FALSE.
1499    IF ( bc_ns == 'dirichlet/radiation' )  bc_ns_dirrad = .TRUE.
1500    IF ( bc_ns == 'radiation/dirichlet' )  bc_ns_raddir = .TRUE.
1501
1502!
1503!-- Non-cyclic lateral boundaries require the multigrid method and Piascek-
1504!-- Willimas or Wicker - Skamarock advection scheme. Several schemes
1505!-- and tools do not work with non-cyclic boundary conditions.
1506    IF ( bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic' )  THEN
1507       IF ( psolver /= 'multigrid' )  THEN
1508          message_string = 'non-cyclic lateral boundaries do not allow ' // &
1509                           'psolver = "' // TRIM( psolver ) // '"'
1510          CALL message( 'check_parameters', 'PA0051', 1, 2, 0, 6, 0 )
1511       ENDIF
1512       IF ( momentum_advec /= 'pw-scheme' .AND. &
1513            momentum_advec /= 'ws-scheme')  THEN
1514          message_string = 'non-cyclic lateral boundaries do not allow ' // &
1515                           'momentum_advec = "' // TRIM( momentum_advec ) // '"'
1516          CALL message( 'check_parameters', 'PA0052', 1, 2, 0, 6, 0 )
1517       ENDIF
1518       IF ( scalar_advec /= 'pw-scheme' .AND. &
1519            scalar_advec /= 'ws-scheme' )  THEN
1520          message_string = 'non-cyclic lateral boundaries do not allow ' // &
1521                           'scalar_advec = "' // TRIM( scalar_advec ) // '"'
1522          CALL message( 'check_parameters', 'PA0053', 1, 2, 0, 6, 0 )
1523       ENDIF
1524       IF ( galilei_transformation )  THEN
1525          message_string = 'non-cyclic lateral boundaries do not allow ' // &
1526                           'galilei_transformation = .T.'
1527          CALL message( 'check_parameters', 'PA0054', 1, 2, 0, 6, 0 )
1528       ENDIF
1529    ENDIF
1530
1531!
1532!-- Bottom boundary condition for the turbulent Kinetic energy
1533    IF ( bc_e_b == 'neumann' )  THEN
1534       ibc_e_b = 1
1535    ELSEIF ( bc_e_b == '(u*)**2+neumann' )  THEN
1536       ibc_e_b = 2
1537       IF ( .NOT. prandtl_layer )  THEN
1538          bc_e_b = 'neumann'
1539          ibc_e_b = 1
1540          message_string = 'boundary condition bc_e_b changed to "' // &
1541                           TRIM( bc_e_b ) // '"'
1542          CALL message( 'check_parameters', 'PA0057', 0, 1, 0, 6, 0 )
1543       ENDIF
1544    ELSE
1545       message_string = 'unknown boundary condition: bc_e_b = "' // &
1546                        TRIM( bc_e_b ) // '"'
1547       CALL message( 'check_parameters', 'PA0058', 1, 2, 0, 6, 0 )
1548    ENDIF
1549
1550!
1551!-- Boundary conditions for perturbation pressure
1552    IF ( bc_p_b == 'dirichlet' )  THEN
1553       ibc_p_b = 0
1554    ELSEIF ( bc_p_b == 'neumann' )  THEN
1555       ibc_p_b = 1
1556    ELSE
1557       message_string = 'unknown boundary condition: bc_p_b = "' // &
1558                        TRIM( bc_p_b ) // '"'
1559       CALL message( 'check_parameters', 'PA0059', 1, 2, 0, 6, 0 )
1560    ENDIF
1561
1562    IF ( bc_p_t == 'dirichlet' )  THEN
1563       ibc_p_t = 0
1564    ELSEIF ( bc_p_t == 'neumann' )  THEN
1565       ibc_p_t = 1
1566    ELSE
1567       message_string = 'unknown boundary condition: bc_p_t = "' // &
1568                        TRIM( bc_p_t ) // '"'
1569       CALL message( 'check_parameters', 'PA0061', 1, 2, 0, 6, 0 )
1570    ENDIF
1571
1572!
1573!-- Boundary conditions for potential temperature
1574    IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
1575       ibc_pt_b = 2
1576    ELSE
1577       IF ( bc_pt_b == 'dirichlet' )  THEN
1578          ibc_pt_b = 0
1579       ELSEIF ( bc_pt_b == 'neumann' )  THEN
1580          ibc_pt_b = 1
1581       ELSE
1582          message_string = 'unknown boundary condition: bc_pt_b = "' // &
1583                           TRIM( bc_pt_b ) // '"'
1584          CALL message( 'check_parameters', 'PA0062', 1, 2, 0, 6, 0 )
1585       ENDIF
1586    ENDIF
1587
1588    IF ( bc_pt_t == 'dirichlet' )  THEN
1589       ibc_pt_t = 0
1590    ELSEIF ( bc_pt_t == 'neumann' )  THEN
1591       ibc_pt_t = 1
1592    ELSEIF ( bc_pt_t == 'initial_gradient' )  THEN
1593       ibc_pt_t = 2
1594    ELSE
1595       message_string = 'unknown boundary condition: bc_pt_t = "' // &
1596                        TRIM( bc_pt_t ) // '"'
1597       CALL message( 'check_parameters', 'PA0063', 1, 2, 0, 6, 0 )
1598    ENDIF
1599
1600    IF ( surface_heatflux == 9999999.9_wp  )  THEN
1601       constant_heatflux = .FALSE.
1602       IF ( large_scale_forcing )  THEN
1603          IF ( ibc_pt_b == 0 )  THEN
1604             constant_heatflux = .FALSE.
1605          ELSEIF ( ibc_pt_b == 1 )  THEN
1606             constant_heatflux = .TRUE.
1607             IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1608                surface_heatflux = shf_surf(1)
1609             ENDIF
1610          ENDIF
1611       ENDIF
1612    ELSE
1613        constant_heatflux = .TRUE.
1614        IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. &
1615               large_scale_forcing ) THEN
1616           surface_heatflux = shf_surf(1)
1617        ENDIF
1618    ENDIF
1619
1620    IF ( top_heatflux     == 9999999.9_wp )  constant_top_heatflux = .FALSE.
1621
1622    IF ( neutral )  THEN
1623
1624       IF ( surface_heatflux /= 0.0_wp  .AND.  surface_heatflux /= 9999999.9_wp ) &
1625       THEN
1626          message_string = 'heatflux must not be set for pure neutral flow'
1627          CALL message( 'check_parameters', 'PA0351', 1, 2, 0, 6, 0 )
1628       ENDIF
1629
1630       IF ( top_heatflux /= 0.0  .AND.  top_heatflux /= 9999999.9_wp ) &
1631       THEN
1632          message_string = 'heatflux must not be set for pure neutral flow'
1633          CALL message( 'check_parameters', 'PA0351', 1, 2, 0, 6, 0 )
1634       ENDIF
1635
1636    ENDIF
1637
1638    IF ( top_momentumflux_u /= 9999999.9_wp  .AND.  &
1639         top_momentumflux_v /= 9999999.9_wp )  THEN
1640       constant_top_momentumflux = .TRUE.
1641    ELSEIF (  .NOT. ( top_momentumflux_u == 9999999.9_wp  .AND.  &
1642           top_momentumflux_v == 9999999.9_wp ) )  THEN
1643       message_string = 'both, top_momentumflux_u AND top_momentumflux_v ' // &
1644                        'must be set'
1645       CALL message( 'check_parameters', 'PA0064', 1, 2, 0, 6, 0 )
1646    ENDIF
1647
1648!
1649!-- A given surface temperature implies Dirichlet boundary condition for
1650!-- temperature. In this case specification of a constant heat flux is
1651!-- forbidden.
1652    IF ( ibc_pt_b == 0  .AND.   constant_heatflux  .AND. &
1653         surface_heatflux /= 0.0_wp )  THEN
1654       message_string = 'boundary_condition: bc_pt_b = "' // TRIM( bc_pt_b ) //&
1655                        '& is not allowed with constant_heatflux = .TRUE.'
1656       CALL message( 'check_parameters', 'PA0065', 1, 2, 0, 6, 0 )
1657    ENDIF
1658    IF ( constant_heatflux  .AND.  pt_surface_initial_change /= 0.0_wp )  THEN
1659       WRITE ( message_string, * )  'constant_heatflux = .TRUE. is not allo', &
1660               'wed with pt_surface_initial_change (/=0) = ', &
1661               pt_surface_initial_change
1662       CALL message( 'check_parameters', 'PA0066', 1, 2, 0, 6, 0 )
1663    ENDIF
1664
1665!
1666!-- A given temperature at the top implies Dirichlet boundary condition for
1667!-- temperature. In this case specification of a constant heat flux is
1668!-- forbidden.
1669    IF ( ibc_pt_t == 0  .AND.   constant_top_heatflux  .AND. &
1670         top_heatflux /= 0.0_wp )  THEN
1671       message_string = 'boundary_condition: bc_pt_t = "' // TRIM( bc_pt_t ) //&
1672                        '" is not allowed with constant_top_heatflux = .TRUE.'
1673       CALL message( 'check_parameters', 'PA0067', 1, 2, 0, 6, 0 )
1674    ENDIF
1675
1676!
1677!-- Boundary conditions for salinity
1678    IF ( ocean )  THEN
1679       IF ( bc_sa_t == 'dirichlet' )  THEN
1680          ibc_sa_t = 0
1681       ELSEIF ( bc_sa_t == 'neumann' )  THEN
1682          ibc_sa_t = 1
1683       ELSE
1684          message_string = 'unknown boundary condition: bc_sa_t = "' // &
1685                           TRIM( bc_sa_t ) // '"'
1686          CALL message( 'check_parameters', 'PA0068', 1, 2, 0, 6, 0 )
1687       ENDIF
1688
1689       IF ( top_salinityflux == 9999999.9_wp )  constant_top_salinityflux = .FALSE.
1690       IF ( ibc_sa_t == 1  .AND.   top_salinityflux == 9999999.9_wp )  THEN
1691          message_string = 'boundary condition: bc_sa_t = "' // &
1692                           TRIM( bc_sa_t ) // '" requires to set ' // &
1693                           'top_salinityflux'
1694          CALL message( 'check_parameters', 'PA0069', 1, 2, 0, 6, 0 )
1695       ENDIF
1696
1697!
1698!--    A fixed salinity at the top implies Dirichlet boundary condition for
1699!--    salinity. In this case specification of a constant salinity flux is
1700!--    forbidden.
1701       IF ( ibc_sa_t == 0  .AND.   constant_top_salinityflux  .AND. &
1702            top_salinityflux /= 0.0_wp )  THEN
1703          message_string = 'boundary condition: bc_sa_t = "' // &
1704                           TRIM( bc_sa_t ) // '" is not allowed with ' // &
1705                           'constant_top_salinityflux = .TRUE.'
1706          CALL message( 'check_parameters', 'PA0070', 1, 2, 0, 6, 0 )
1707       ENDIF
1708
1709    ENDIF
1710
1711!
1712!-- In case of humidity or passive scalar, set boundary conditions for total
1713!-- water content / scalar
1714    IF ( humidity  .OR.  passive_scalar ) THEN
1715       IF ( humidity )  THEN
1716          sq = 'q'
1717       ELSE
1718          sq = 's'
1719       ENDIF
1720       IF ( bc_q_b == 'dirichlet' )  THEN
1721          ibc_q_b = 0
1722       ELSEIF ( bc_q_b == 'neumann' )  THEN
1723          ibc_q_b = 1
1724       ELSE
1725          message_string = 'unknown boundary condition: bc_' // TRIM( sq ) // &
1726                           '_b ="' // TRIM( bc_q_b ) // '"'
1727          CALL message( 'check_parameters', 'PA0071', 1, 2, 0, 6, 0 )
1728       ENDIF
1729       IF ( bc_q_t == 'dirichlet' )  THEN
1730          ibc_q_t = 0
1731       ELSEIF ( bc_q_t == 'neumann' )  THEN
1732          ibc_q_t = 1
1733       ELSE
1734          message_string = 'unknown boundary condition: bc_' // TRIM( sq ) // &
1735                           '_t ="' // TRIM( bc_q_t ) // '"'
1736          CALL message( 'check_parameters', 'PA0072', 1, 2, 0, 6, 0 )
1737       ENDIF
1738
1739       IF ( surface_waterflux == 9999999.9_wp  )  THEN
1740          constant_waterflux = .FALSE.
1741          IF ( large_scale_forcing )  THEN
1742             IF ( ibc_q_b == 0 )  THEN
1743                constant_waterflux = .FALSE.
1744             ELSEIF ( ibc_q_b == 1 )  THEN
1745                constant_waterflux = .TRUE.
1746                IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1747                   surface_waterflux = qsws_surf(1)
1748                ENDIF
1749             ENDIF
1750          ENDIF
1751       ELSE
1752          constant_waterflux = .TRUE.
1753          IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. &
1754                 large_scale_forcing ) THEN
1755             surface_waterflux = qsws_surf(1)
1756          ENDIF
1757       ENDIF
1758
1759!
1760!--    A given surface humidity implies Dirichlet boundary condition for
1761!--    humidity. In this case specification of a constant water flux is
1762!--    forbidden.
1763       IF ( ibc_q_b == 0  .AND.  constant_waterflux )  THEN
1764          message_string = 'boundary condition: bc_' // TRIM( sq ) // '_b ' // &
1765                           '= "' // TRIM( bc_q_b ) // '" is not allowed wi' // &
1766                           'th prescribed surface flux'
1767          CALL message( 'check_parameters', 'PA0073', 1, 2, 0, 6, 0 )
1768       ENDIF
1769       IF ( constant_waterflux  .AND.  q_surface_initial_change /= 0.0_wp )  THEN
1770          WRITE( message_string, * )  'a prescribed surface flux is not allo', &
1771                 'wed with ', sq, '_surface_initial_change (/=0) = ', &
1772                 q_surface_initial_change
1773          CALL message( 'check_parameters', 'PA0074', 1, 2, 0, 6, 0 )
1774       ENDIF
1775
1776    ENDIF
1777!
1778!-- Boundary conditions for horizontal components of wind speed
1779    IF ( bc_uv_b == 'dirichlet' )  THEN
1780       ibc_uv_b = 0
1781    ELSEIF ( bc_uv_b == 'neumann' )  THEN
1782       ibc_uv_b = 1
1783       IF ( prandtl_layer )  THEN
1784          message_string = 'boundary condition: bc_uv_b = "' // &
1785               TRIM( bc_uv_b ) // '" is not allowed with prandtl_layer = .TRUE.'
1786          CALL message( 'check_parameters', 'PA0075', 1, 2, 0, 6, 0 )
1787       ENDIF
1788    ELSE
1789       message_string = 'unknown boundary condition: bc_uv_b = "' // &
1790                        TRIM( bc_uv_b ) // '"'
1791       CALL message( 'check_parameters', 'PA0076', 1, 2, 0, 6, 0 )
1792    ENDIF
1793!
1794!-- In case of coupled simulations u and v at the ground in atmosphere will be
1795!-- assigned with the u and v values of the ocean surface
1796    IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
1797       ibc_uv_b = 2
1798    ENDIF
1799
1800    IF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
1801       bc_uv_t = 'neumann'
1802       ibc_uv_t = 1
1803    ELSE
1804       IF ( bc_uv_t == 'dirichlet' .OR. bc_uv_t == 'dirichlet_0' )  THEN
1805          ibc_uv_t = 0
1806          IF ( bc_uv_t == 'dirichlet_0' )  THEN
1807!
1808!--          Velocities for the initial u,v-profiles are set zero at the top
1809!--          in case of dirichlet_0 conditions
1810             u_init(nzt+1)    = 0.0_wp
1811             v_init(nzt+1)    = 0.0_wp
1812          ENDIF
1813       ELSEIF ( bc_uv_t == 'neumann' )  THEN
1814          ibc_uv_t = 1
1815       ELSE
1816          message_string = 'unknown boundary condition: bc_uv_t = "' // &
1817                           TRIM( bc_uv_t ) // '"'
1818          CALL message( 'check_parameters', 'PA0077', 1, 2, 0, 6, 0 )
1819       ENDIF
1820    ENDIF
1821
1822!
1823!-- Compute and check, respectively, the Rayleigh Damping parameter
1824    IF ( rayleigh_damping_factor == -1.0_wp )  THEN
1825       rayleigh_damping_factor = 0.0_wp
1826    ELSE
1827       IF ( rayleigh_damping_factor < 0.0 .OR. rayleigh_damping_factor > 1.0_wp ) &
1828       THEN
1829          WRITE( message_string, * )  'rayleigh_damping_factor = ', &
1830                              rayleigh_damping_factor, ' out of range [0.0,1.0]'
1831          CALL message( 'check_parameters', 'PA0078', 1, 2, 0, 6, 0 )
1832       ENDIF
1833    ENDIF
1834
1835    IF ( rayleigh_damping_height == -1.0_wp )  THEN
1836       IF ( .NOT. ocean )  THEN
1837          rayleigh_damping_height = 0.66666666666_wp * zu(nzt)
1838       ELSE
1839          rayleigh_damping_height = 0.66666666666_wp * zu(nzb)
1840       ENDIF
1841    ELSE
1842       IF ( .NOT. ocean )  THEN
1843          IF ( rayleigh_damping_height < 0.0_wp  .OR. &
1844               rayleigh_damping_height > zu(nzt) )  THEN
1845             WRITE( message_string, * )  'rayleigh_damping_height = ', &
1846                   rayleigh_damping_height, ' out of range [0.0,', zu(nzt), ']'
1847             CALL message( 'check_parameters', 'PA0079', 1, 2, 0, 6, 0 )
1848          ENDIF
1849       ELSE
1850          IF ( rayleigh_damping_height > 0.0_wp  .OR. &
1851               rayleigh_damping_height < zu(nzb) )  THEN
1852             WRITE( message_string, * )  'rayleigh_damping_height = ', &
1853                   rayleigh_damping_height, ' out of range [0.0,', zu(nzb), ']'
1854             CALL message( 'check_parameters', 'PA0079', 1, 2, 0, 6, 0 )
1855          ENDIF
1856       ENDIF
1857    ENDIF
1858
1859!
1860!-- Check number of chosen statistic regions. More than 10 regions are not
1861!-- allowed, because so far no more than 10 corresponding output files can
1862!-- be opened (cf. check_open)
1863    IF ( statistic_regions > 9  .OR.  statistic_regions < 0 )  THEN
1864       WRITE ( message_string, * ) 'number of statistic_regions = ', &
1865                   statistic_regions+1, ' but only 10 regions are allowed'
1866       CALL message( 'check_parameters', 'PA0082', 1, 2, 0, 6, 0 )
1867    ENDIF
1868    IF ( normalizing_region > statistic_regions  .OR. &
1869         normalizing_region < 0)  THEN
1870       WRITE ( message_string, * ) 'normalizing_region = ', &
1871                normalizing_region, ' must be >= 0 and <= ',statistic_regions, &
1872                ' (value of statistic_regions)'
1873       CALL message( 'check_parameters', 'PA0083', 1, 2, 0, 6, 0 )
1874    ENDIF
1875
1876!
1877!-- Set the default intervals for data output, if necessary
1878!-- NOTE: dt_dosp has already been set in package_parin
1879    IF ( dt_data_output /= 9999999.9_wp )  THEN
1880       IF ( dt_dopr           == 9999999.9_wp )  dt_dopr           = dt_data_output
1881       IF ( dt_dopts          == 9999999.9_wp )  dt_dopts          = dt_data_output
1882       IF ( dt_do2d_xy        == 9999999.9_wp )  dt_do2d_xy        = dt_data_output
1883       IF ( dt_do2d_xz        == 9999999.9_wp )  dt_do2d_xz        = dt_data_output
1884       IF ( dt_do2d_yz        == 9999999.9_wp )  dt_do2d_yz        = dt_data_output
1885       IF ( dt_do3d           == 9999999.9_wp )  dt_do3d           = dt_data_output
1886       IF ( dt_data_output_av == 9999999.9_wp )  dt_data_output_av = dt_data_output
1887       DO  mid = 1, max_masks
1888          IF ( dt_domask(mid) == 9999999.9_wp )  dt_domask(mid)    = dt_data_output
1889       ENDDO
1890    ENDIF
1891
1892!
1893!-- Set the default skip time intervals for data output, if necessary
1894    IF ( skip_time_dopr    == 9999999.9_wp ) &
1895                                       skip_time_dopr    = skip_time_data_output
1896    IF ( skip_time_dosp    == 9999999.9_wp ) &
1897                                       skip_time_dosp    = skip_time_data_output
1898    IF ( skip_time_do2d_xy == 9999999.9_wp ) &
1899                                       skip_time_do2d_xy = skip_time_data_output
1900    IF ( skip_time_do2d_xz == 9999999.9_wp ) &
1901                                       skip_time_do2d_xz = skip_time_data_output
1902    IF ( skip_time_do2d_yz == 9999999.9_wp ) &
1903                                       skip_time_do2d_yz = skip_time_data_output
1904    IF ( skip_time_do3d    == 9999999.9_wp ) &
1905                                       skip_time_do3d    = skip_time_data_output
1906    IF ( skip_time_data_output_av == 9999999.9_wp ) &
1907                                skip_time_data_output_av = skip_time_data_output
1908    DO  mid = 1, max_masks
1909       IF ( skip_time_domask(mid) == 9999999.9_wp ) &
1910                                skip_time_domask(mid)    = skip_time_data_output
1911    ENDDO
1912
1913!
1914!-- Check the average intervals (first for 3d-data, then for profiles and
1915!-- spectra)
1916    IF ( averaging_interval > dt_data_output_av )  THEN
1917       WRITE( message_string, * )  'averaging_interval = ', &
1918             averaging_interval, ' must be <= dt_data_output = ', dt_data_output
1919       CALL message( 'check_parameters', 'PA0085', 1, 2, 0, 6, 0 )
1920    ENDIF
1921
1922    IF ( averaging_interval_pr == 9999999.9_wp )  THEN
1923       averaging_interval_pr = averaging_interval
1924    ENDIF
1925
1926    IF ( averaging_interval_pr > dt_dopr )  THEN
1927       WRITE( message_string, * )  'averaging_interval_pr = ', &
1928             averaging_interval_pr, ' must be <= dt_dopr = ', dt_dopr
1929       CALL message( 'check_parameters', 'PA0086', 1, 2, 0, 6, 0 )
1930    ENDIF
1931
1932    IF ( averaging_interval_sp == 9999999.9_wp )  THEN
1933       averaging_interval_sp = averaging_interval
1934    ENDIF
1935
1936    IF ( averaging_interval_sp > dt_dosp )  THEN
1937       WRITE( message_string, * )  'averaging_interval_sp = ', &
1938             averaging_interval_sp, ' must be <= dt_dosp = ', dt_dosp
1939       CALL message( 'check_parameters', 'PA0087', 1, 2, 0, 6, 0 )
1940    ENDIF
1941
1942!
1943!-- Set the default interval for profiles entering the temporal average
1944    IF ( dt_averaging_input_pr == 9999999.9_wp )  THEN
1945       dt_averaging_input_pr = dt_averaging_input
1946    ENDIF
1947
1948!
1949!-- Set the default interval for the output of timeseries to a reasonable
1950!-- value (tries to minimize the number of calls of flow_statistics)
1951    IF ( dt_dots == 9999999.9_wp )  THEN
1952       IF ( averaging_interval_pr == 0.0_wp )  THEN
1953          dt_dots = MIN( dt_run_control, dt_dopr )
1954       ELSE
1955          dt_dots = MIN( dt_run_control, dt_averaging_input_pr )
1956       ENDIF
1957    ENDIF
1958
1959!
1960!-- Check the sample rate for averaging (first for 3d-data, then for profiles)
1961    IF ( dt_averaging_input > averaging_interval )  THEN
1962       WRITE( message_string, * )  'dt_averaging_input = ', &
1963                dt_averaging_input, ' must be <= averaging_interval = ', &
1964                averaging_interval
1965       CALL message( 'check_parameters', 'PA0088', 1, 2, 0, 6, 0 )
1966    ENDIF
1967
1968    IF ( dt_averaging_input_pr > averaging_interval_pr )  THEN
1969       WRITE( message_string, * )  'dt_averaging_input_pr = ', &
1970                dt_averaging_input_pr, ' must be <= averaging_interval_pr = ', &
1971                averaging_interval_pr
1972       CALL message( 'check_parameters', 'PA0089', 1, 2, 0, 6, 0 )
1973    ENDIF
1974
1975!
1976!-- Set the default value for the integration interval of precipitation amount
1977    IF ( precipitation )  THEN
1978       IF ( precipitation_amount_interval == 9999999.9_wp )  THEN
1979          precipitation_amount_interval = dt_do2d_xy
1980       ELSE
1981          IF ( precipitation_amount_interval > dt_do2d_xy )  THEN
1982             WRITE( message_string, * )  'precipitation_amount_interval = ', &
1983                 precipitation_amount_interval, ' must not be larger than ', &
1984                 'dt_do2d_xy = ', dt_do2d_xy
1985             CALL message( 'check_parameters', 'PA0090', 1, 2, 0, 6, 0 )
1986          ENDIF
1987       ENDIF
1988    ENDIF
1989
1990!
1991!-- Determine the number of output profiles and check whether they are
1992!-- permissible
1993    DO  WHILE ( data_output_pr(dopr_n+1) /= '          ' )
1994
1995       dopr_n = dopr_n + 1
1996       i = dopr_n
1997
1998!
1999!--    Determine internal profile number (for hom, homs)
2000!--    and store height levels
2001       SELECT CASE ( TRIM( data_output_pr(i) ) )
2002
2003          CASE ( 'u', '#u' )
2004             dopr_index(i) = 1
2005             dopr_unit(i)  = 'm/s'
2006             hom(:,2,1,:)  = SPREAD( zu, 2, statistic_regions+1 )
2007             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2008                dopr_initial_index(i) = 5
2009                hom(:,2,5,:)          = SPREAD( zu, 2, statistic_regions+1 )
2010                data_output_pr(i)     = data_output_pr(i)(2:)
2011             ENDIF
2012
2013          CASE ( 'v', '#v' )
2014             dopr_index(i) = 2
2015             dopr_unit(i)  = 'm/s'
2016             hom(:,2,2,:)  = SPREAD( zu, 2, statistic_regions+1 )
2017             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2018                dopr_initial_index(i) = 6
2019                hom(:,2,6,:)          = SPREAD( zu, 2, statistic_regions+1 )
2020                data_output_pr(i)     = data_output_pr(i)(2:)
2021             ENDIF
2022
2023          CASE ( 'w' )
2024             dopr_index(i) = 3
2025             dopr_unit(i)  = 'm/s'
2026             hom(:,2,3,:)  = SPREAD( zw, 2, statistic_regions+1 )
2027
2028          CASE ( 'pt', '#pt' )
2029             IF ( .NOT. cloud_physics ) THEN
2030                dopr_index(i) = 4
2031                dopr_unit(i)  = 'K'
2032                hom(:,2,4,:)  = SPREAD( zu, 2, statistic_regions+1 )
2033                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2034                   dopr_initial_index(i) = 7
2035                   hom(:,2,7,:)          = SPREAD( zu, 2, statistic_regions+1 )
2036                   hom(nzb,2,7,:)        = 0.0_wp    ! because zu(nzb) is negative
2037                   data_output_pr(i)     = data_output_pr(i)(2:)
2038                ENDIF
2039             ELSE
2040                dopr_index(i) = 43
2041                dopr_unit(i)  = 'K'
2042                hom(:,2,43,:)  = SPREAD( zu, 2, statistic_regions+1 )
2043                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2044                   dopr_initial_index(i) = 28
2045                   hom(:,2,28,:)         = SPREAD( zu, 2, statistic_regions+1 )
2046                   hom(nzb,2,28,:)       = 0.0_wp    ! because zu(nzb) is negative
2047                   data_output_pr(i)     = data_output_pr(i)(2:)
2048                ENDIF
2049             ENDIF
2050
2051          CASE ( 'e' )
2052             dopr_index(i)  = 8
2053             dopr_unit(i)   = 'm2/s2'
2054             hom(:,2,8,:)   = SPREAD( zu, 2, statistic_regions+1 )
2055             hom(nzb,2,8,:) = 0.0_wp
2056
2057          CASE ( 'km', '#km' )
2058             dopr_index(i)  = 9
2059             dopr_unit(i)   = 'm2/s'
2060             hom(:,2,9,:)   = SPREAD( zu, 2, statistic_regions+1 )
2061             hom(nzb,2,9,:) = 0.0_wp
2062             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2063                dopr_initial_index(i) = 23
2064                hom(:,2,23,:)         = hom(:,2,9,:)
2065                data_output_pr(i)     = data_output_pr(i)(2:)
2066             ENDIF
2067
2068          CASE ( 'kh', '#kh' )
2069             dopr_index(i)   = 10
2070             dopr_unit(i)    = 'm2/s'
2071             hom(:,2,10,:)   = SPREAD( zu, 2, statistic_regions+1 )
2072             hom(nzb,2,10,:) = 0.0_wp
2073             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2074                dopr_initial_index(i) = 24
2075                hom(:,2,24,:)         = hom(:,2,10,:)
2076                data_output_pr(i)     = data_output_pr(i)(2:)
2077             ENDIF
2078
2079          CASE ( 'l', '#l' )
2080             dopr_index(i)   = 11
2081             dopr_unit(i)    = 'm'
2082             hom(:,2,11,:)   = SPREAD( zu, 2, statistic_regions+1 )
2083             hom(nzb,2,11,:) = 0.0_wp
2084             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2085                dopr_initial_index(i) = 25
2086                hom(:,2,25,:)         = hom(:,2,11,:)
2087                data_output_pr(i)     = data_output_pr(i)(2:)
2088             ENDIF
2089
2090          CASE ( 'w"u"' )
2091             dopr_index(i) = 12
2092             dopr_unit(i)  = 'm2/s2'
2093             hom(:,2,12,:) = SPREAD( zw, 2, statistic_regions+1 )
2094             IF ( prandtl_layer )  hom(nzb,2,12,:) = zu(1)
2095
2096          CASE ( 'w*u*' )
2097             dopr_index(i) = 13
2098             dopr_unit(i)  = 'm2/s2'
2099             hom(:,2,13,:) = SPREAD( zw, 2, statistic_regions+1 )
2100
2101          CASE ( 'w"v"' )
2102             dopr_index(i) = 14
2103             dopr_unit(i)  = 'm2/s2'
2104             hom(:,2,14,:) = SPREAD( zw, 2, statistic_regions+1 )
2105             IF ( prandtl_layer )  hom(nzb,2,14,:) = zu(1)
2106
2107          CASE ( 'w*v*' )
2108             dopr_index(i) = 15
2109             dopr_unit(i)  = 'm2/s2'
2110             hom(:,2,15,:) = SPREAD( zw, 2, statistic_regions+1 )
2111
2112          CASE ( 'w"pt"' )
2113             dopr_index(i) = 16
2114             dopr_unit(i)  = 'K m/s'
2115             hom(:,2,16,:) = SPREAD( zw, 2, statistic_regions+1 )
2116
2117          CASE ( 'w*pt*' )
2118             dopr_index(i) = 17
2119             dopr_unit(i)  = 'K m/s'
2120             hom(:,2,17,:) = SPREAD( zw, 2, statistic_regions+1 )
2121
2122          CASE ( 'wpt' )
2123             dopr_index(i) = 18
2124             dopr_unit(i)  = 'K m/s'
2125             hom(:,2,18,:) = SPREAD( zw, 2, statistic_regions+1 )
2126
2127          CASE ( 'wu' )
2128             dopr_index(i) = 19
2129             dopr_unit(i)  = 'm2/s2'
2130             hom(:,2,19,:) = SPREAD( zw, 2, statistic_regions+1 )
2131             IF ( prandtl_layer )  hom(nzb,2,19,:) = zu(1)
2132
2133          CASE ( 'wv' )
2134             dopr_index(i) = 20
2135             dopr_unit(i)  = 'm2/s2'
2136             hom(:,2,20,:) = SPREAD( zw, 2, statistic_regions+1 )
2137             IF ( prandtl_layer )  hom(nzb,2,20,:) = zu(1)
2138
2139          CASE ( 'w*pt*BC' )
2140             dopr_index(i) = 21
2141             dopr_unit(i)  = 'K m/s'
2142             hom(:,2,21,:) = SPREAD( zw, 2, statistic_regions+1 )
2143
2144          CASE ( 'wptBC' )
2145             dopr_index(i) = 22
2146             dopr_unit(i)  = 'K m/s'
2147             hom(:,2,22,:) = SPREAD( zw, 2, statistic_regions+1 )
2148
2149          CASE ( 'sa', '#sa' )
2150             IF ( .NOT. ocean )  THEN
2151                message_string = 'data_output_pr = ' // &
2152                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2153                                 'lemented for ocean = .FALSE.'
2154                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2155             ELSE
2156                dopr_index(i) = 23
2157                dopr_unit(i)  = 'psu'
2158                hom(:,2,23,:) = SPREAD( zu, 2, statistic_regions+1 )
2159                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2160                   dopr_initial_index(i) = 26
2161                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
2162                   hom(nzb,2,26,:)       = 0.0_wp    ! because zu(nzb) is negative
2163                   data_output_pr(i)     = data_output_pr(i)(2:)
2164                ENDIF
2165             ENDIF
2166
2167          CASE ( 'u*2' )
2168             dopr_index(i) = 30
2169             dopr_unit(i)  = 'm2/s2'
2170             hom(:,2,30,:) = SPREAD( zu, 2, statistic_regions+1 )
2171
2172          CASE ( 'v*2' )
2173             dopr_index(i) = 31
2174             dopr_unit(i)  = 'm2/s2'
2175             hom(:,2,31,:) = SPREAD( zu, 2, statistic_regions+1 )
2176
2177          CASE ( 'w*2' )
2178             dopr_index(i) = 32
2179             dopr_unit(i)  = 'm2/s2'
2180             hom(:,2,32,:) = SPREAD( zw, 2, statistic_regions+1 )
2181
2182          CASE ( 'pt*2' )
2183             dopr_index(i) = 33
2184             dopr_unit(i)  = 'K2'
2185             hom(:,2,33,:) = SPREAD( zu, 2, statistic_regions+1 )
2186
2187          CASE ( 'e*' )
2188             dopr_index(i) = 34
2189             dopr_unit(i)  = 'm2/s2'
2190             hom(:,2,34,:) = SPREAD( zu, 2, statistic_regions+1 )
2191
2192          CASE ( 'w*2pt*' )
2193             dopr_index(i) = 35
2194             dopr_unit(i)  = 'K m2/s2'
2195             hom(:,2,35,:) = SPREAD( zw, 2, statistic_regions+1 )
2196
2197          CASE ( 'w*pt*2' )
2198             dopr_index(i) = 36
2199             dopr_unit(i)  = 'K2 m/s'
2200             hom(:,2,36,:) = SPREAD( zw, 2, statistic_regions+1 )
2201
2202          CASE ( 'w*e*' )
2203             dopr_index(i) = 37
2204             dopr_unit(i)  = 'm3/s3'
2205             hom(:,2,37,:) = SPREAD( zw, 2, statistic_regions+1 )
2206
2207          CASE ( 'w*3' )
2208             dopr_index(i) = 38
2209             dopr_unit(i)  = 'm3/s3'
2210             hom(:,2,38,:) = SPREAD( zw, 2, statistic_regions+1 )
2211
2212          CASE ( 'Sw' )
2213             dopr_index(i) = 39
2214             dopr_unit(i)  = 'none'
2215             hom(:,2,39,:) = SPREAD( zw, 2, statistic_regions+1 )
2216
2217          CASE ( 'p' )
2218             dopr_index(i) = 40
2219             dopr_unit(i)  = 'Pa'
2220             hom(:,2,40,:) = SPREAD( zu, 2, statistic_regions+1 )
2221
2222          CASE ( 'q', '#q' )
2223             IF ( .NOT. humidity )  THEN
2224                message_string = 'data_output_pr = ' // &
2225                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2226                                 'lemented for humidity = .FALSE.'
2227                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2228             ELSE
2229                dopr_index(i) = 41
2230                dopr_unit(i)  = 'kg/kg'
2231                hom(:,2,41,:) = SPREAD( zu, 2, statistic_regions+1 )
2232                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2233                   dopr_initial_index(i) = 26
2234                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
2235                   hom(nzb,2,26,:)       = 0.0_wp    ! because zu(nzb) is negative
2236                   data_output_pr(i)     = data_output_pr(i)(2:)
2237                ENDIF
2238             ENDIF
2239
2240          CASE ( 's', '#s' )
2241             IF ( .NOT. passive_scalar )  THEN
2242                message_string = 'data_output_pr = ' // &
2243                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2244                                 'lemented for passive_scalar = .FALSE.'
2245                CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
2246             ELSE
2247                dopr_index(i) = 41
2248                dopr_unit(i)  = 'kg/m3'
2249                hom(:,2,41,:) = SPREAD( zu, 2, statistic_regions+1 )
2250                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2251                   dopr_initial_index(i) = 26
2252                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
2253                   hom(nzb,2,26,:)       = 0.0_wp    ! because zu(nzb) is negative
2254                   data_output_pr(i)     = data_output_pr(i)(2:)
2255                ENDIF
2256             ENDIF
2257
2258          CASE ( 'qv', '#qv' )
2259             IF ( .NOT. cloud_physics ) THEN
2260                dopr_index(i) = 41
2261                dopr_unit(i)  = 'kg/kg'
2262                hom(:,2,41,:) = SPREAD( zu, 2, statistic_regions+1 )
2263                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2264                   dopr_initial_index(i) = 26
2265                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
2266                   hom(nzb,2,26,:)       = 0.0_wp    ! because zu(nzb) is negative
2267                   data_output_pr(i)     = data_output_pr(i)(2:)
2268                ENDIF
2269             ELSE
2270                dopr_index(i) = 42
2271                dopr_unit(i)  = 'kg/kg'
2272                hom(:,2,42,:) = SPREAD( zu, 2, statistic_regions+1 )
2273                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2274                   dopr_initial_index(i) = 27
2275                   hom(:,2,27,:)         = SPREAD( zu, 2, statistic_regions+1 )
2276                   hom(nzb,2,27,:)       = 0.0_wp   ! because zu(nzb) is negative
2277                   data_output_pr(i)     = data_output_pr(i)(2:)
2278                ENDIF
2279             ENDIF
2280
2281          CASE ( 'lpt', '#lpt' )
2282             IF ( .NOT. cloud_physics ) THEN
2283                message_string = 'data_output_pr = ' // &
2284                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2285                                 'lemented for cloud_physics = .FALSE.'
2286                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
2287             ELSE
2288                dopr_index(i) = 4
2289                dopr_unit(i)  = 'K'
2290                hom(:,2,4,:)  = SPREAD( zu, 2, statistic_regions+1 )
2291                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2292                   dopr_initial_index(i) = 7
2293                   hom(:,2,7,:)          = SPREAD( zu, 2, statistic_regions+1 )
2294                   hom(nzb,2,7,:)        = 0.0_wp    ! because zu(nzb) is negative
2295                   data_output_pr(i)     = data_output_pr(i)(2:)
2296                ENDIF
2297             ENDIF
2298
2299          CASE ( 'vpt', '#vpt' )
2300             dopr_index(i) = 44
2301             dopr_unit(i)  = 'K'
2302             hom(:,2,44,:) = SPREAD( zu, 2, statistic_regions+1 )
2303             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2304                dopr_initial_index(i) = 29
2305                hom(:,2,29,:)         = SPREAD( zu, 2, statistic_regions+1 )
2306                hom(nzb,2,29,:)       = 0.0_wp    ! because zu(nzb) is negative
2307                data_output_pr(i)     = data_output_pr(i)(2:)
2308             ENDIF
2309
2310          CASE ( 'w"vpt"' )
2311             dopr_index(i) = 45
2312             dopr_unit(i)  = 'K m/s'
2313             hom(:,2,45,:) = SPREAD( zw, 2, statistic_regions+1 )
2314
2315          CASE ( 'w*vpt*' )
2316             dopr_index(i) = 46
2317             dopr_unit(i)  = 'K m/s'
2318             hom(:,2,46,:) = SPREAD( zw, 2, statistic_regions+1 )
2319
2320          CASE ( 'wvpt' )
2321             dopr_index(i) = 47
2322             dopr_unit(i)  = 'K m/s'
2323             hom(:,2,47,:) = SPREAD( zw, 2, statistic_regions+1 )
2324
2325          CASE ( 'w"q"' )
2326             IF ( .NOT. humidity )  THEN
2327                message_string = 'data_output_pr = ' // &
2328                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2329                                 'lemented for humidity = .FALSE.'
2330                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2331             ELSE
2332                dopr_index(i) = 48
2333                dopr_unit(i)  = 'kg/kg m/s'
2334                hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
2335             ENDIF
2336
2337          CASE ( 'w*q*' )
2338             IF ( .NOT. humidity )  THEN
2339                message_string = 'data_output_pr = ' // &
2340                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2341                                 'lemented for humidity = .FALSE.'
2342                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2343             ELSE
2344                dopr_index(i) = 49
2345                dopr_unit(i)  = 'kg/kg m/s'
2346                hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
2347             ENDIF
2348
2349          CASE ( 'wq' )
2350             IF ( .NOT. humidity )  THEN
2351                message_string = 'data_output_pr = ' // &
2352                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2353                                 'lemented for humidity = .FALSE.'
2354                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2355             ELSE
2356                dopr_index(i) = 50
2357                dopr_unit(i)  = 'kg/kg m/s'
2358                hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
2359             ENDIF
2360
2361          CASE ( 'w"s"' )
2362             IF ( .NOT. passive_scalar ) THEN
2363                message_string = 'data_output_pr = ' // &
2364                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2365                                 'lemented for passive_scalar = .FALSE.'
2366                CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
2367             ELSE
2368                dopr_index(i) = 48
2369                dopr_unit(i)  = 'kg/m3 m/s'
2370                hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
2371             ENDIF
2372
2373          CASE ( 'w*s*' )
2374             IF ( .NOT. passive_scalar ) THEN
2375                message_string = 'data_output_pr = ' // &
2376                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2377                                 'lemented for passive_scalar = .FALSE.'
2378                CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
2379             ELSE
2380                dopr_index(i) = 49
2381                dopr_unit(i)  = 'kg/m3 m/s'
2382                hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
2383             ENDIF
2384
2385          CASE ( 'ws' )
2386             IF ( .NOT. passive_scalar ) THEN
2387                message_string = 'data_output_pr = ' // &
2388                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2389                                 'lemented for passive_scalar = .FALSE.'
2390                CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
2391             ELSE
2392                dopr_index(i) = 50
2393                dopr_unit(i)  = 'kg/m3 m/s'
2394                hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
2395             ENDIF
2396
2397          CASE ( 'w"qv"' )
2398             IF ( humidity  .AND.  .NOT. cloud_physics ) &
2399             THEN
2400                dopr_index(i) = 48
2401                dopr_unit(i)  = 'kg/kg m/s'
2402                hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
2403             ELSEIF( humidity .AND. cloud_physics ) THEN
2404                dopr_index(i) = 51
2405                dopr_unit(i)  = 'kg/kg m/s'
2406                hom(:,2,51,:) = SPREAD( zw, 2, statistic_regions+1 )
2407             ELSE
2408                message_string = 'data_output_pr = ' // &
2409                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2410                                 'lemented for cloud_physics = .FALSE. an&' // &
2411                                 'd humidity = .FALSE.'
2412                CALL message( 'check_parameters', 'PA0095', 1, 2, 0, 6, 0 )
2413             ENDIF
2414
2415          CASE ( 'w*qv*' )
2416             IF ( humidity  .AND.  .NOT. cloud_physics ) &
2417             THEN
2418                dopr_index(i) = 49
2419                dopr_unit(i)  = 'kg/kg m/s'
2420                hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
2421             ELSEIF( humidity .AND. cloud_physics ) THEN
2422                dopr_index(i) = 52
2423                dopr_unit(i)  = 'kg/kg m/s'
2424                hom(:,2,52,:) = SPREAD( zw, 2, statistic_regions+1 )
2425             ELSE
2426                message_string = 'data_output_pr = ' // &
2427                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2428                                 'lemented for cloud_physics = .FALSE. an&' // &
2429                                 'd humidity = .FALSE.'
2430                CALL message( 'check_parameters', 'PA0095', 1, 2, 0, 6, 0 )
2431             ENDIF
2432
2433          CASE ( 'wqv' )
2434             IF ( humidity  .AND.  .NOT. cloud_physics ) &
2435             THEN
2436                dopr_index(i) = 50
2437                dopr_unit(i)  = 'kg/kg m/s'
2438                hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
2439             ELSEIF( humidity .AND. cloud_physics ) THEN
2440                dopr_index(i) = 53
2441                dopr_unit(i)  = 'kg/kg m/s'
2442                hom(:,2,53,:) = SPREAD( zw, 2, statistic_regions+1 )
2443             ELSE
2444                message_string = 'data_output_pr = ' //                        &
2445                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2446                                 'lemented for cloud_physics = .FALSE. an&' // &
2447                                 'd humidity = .FALSE.'
2448                CALL message( 'check_parameters', 'PA0095', 1, 2, 0, 6, 0 )
2449             ENDIF
2450
2451          CASE ( 'ql' )
2452             IF ( .NOT. cloud_physics  .AND.  .NOT. cloud_droplets )  THEN
2453                message_string = 'data_output_pr = ' //                        &
2454                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2455                                 'lemented for cloud_physics = .FALSE. or'  // &
2456                                 '&cloud_droplets = .FALSE.'
2457                CALL message( 'check_parameters', 'PA0096', 1, 2, 0, 6, 0 )
2458             ELSE
2459                dopr_index(i) = 54
2460                dopr_unit(i)  = 'kg/kg'
2461                hom(:,2,54,:)  = SPREAD( zu, 2, statistic_regions+1 )
2462             ENDIF
2463
2464          CASE ( 'w*u*u*:dz' )
2465             dopr_index(i) = 55
2466             dopr_unit(i)  = 'm2/s3'
2467             hom(:,2,55,:) = SPREAD( zu, 2, statistic_regions+1 )
2468
2469          CASE ( 'w*p*:dz' )
2470             dopr_index(i) = 56
2471             dopr_unit(i)  = 'm2/s3'
2472             hom(:,2,56,:) = SPREAD( zw, 2, statistic_regions+1 )
2473
2474          CASE ( 'w"e:dz' )
2475             dopr_index(i) = 57
2476             dopr_unit(i)  = 'm2/s3'
2477             hom(:,2,57,:) = SPREAD( zu, 2, statistic_regions+1 )
2478
2479
2480          CASE ( 'u"pt"' )
2481             dopr_index(i) = 58
2482             dopr_unit(i)  = 'K m/s'
2483             hom(:,2,58,:) = SPREAD( zu, 2, statistic_regions+1 )
2484
2485          CASE ( 'u*pt*' )
2486             dopr_index(i) = 59
2487             dopr_unit(i)  = 'K m/s'
2488             hom(:,2,59,:) = SPREAD( zu, 2, statistic_regions+1 )
2489
2490          CASE ( 'upt_t' )
2491             dopr_index(i) = 60
2492             dopr_unit(i)  = 'K m/s'
2493             hom(:,2,60,:) = SPREAD( zu, 2, statistic_regions+1 )
2494
2495          CASE ( 'v"pt"' )
2496             dopr_index(i) = 61
2497             dopr_unit(i)  = 'K m/s'
2498             hom(:,2,61,:) = SPREAD( zu, 2, statistic_regions+1 )
2499             
2500          CASE ( 'v*pt*' )
2501             dopr_index(i) = 62
2502             dopr_unit(i)  = 'K m/s'
2503             hom(:,2,62,:) = SPREAD( zu, 2, statistic_regions+1 )
2504
2505          CASE ( 'vpt_t' )
2506             dopr_index(i) = 63
2507             dopr_unit(i)  = 'K m/s'
2508             hom(:,2,63,:) = SPREAD( zu, 2, statistic_regions+1 )
2509
2510          CASE ( 'rho' )
2511             IF ( .NOT. ocean ) THEN
2512                message_string = 'data_output_pr = ' //                        &
2513                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2514                                 'lemented for ocean = .FALSE.'
2515                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2516             ELSE
2517                dopr_index(i) = 64
2518                dopr_unit(i)  = 'kg/m3'
2519                hom(:,2,64,:) = SPREAD( zu, 2, statistic_regions+1 )
2520                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2521                   dopr_initial_index(i) = 77
2522                   hom(:,2,77,:)         = SPREAD( zu, 2, statistic_regions+1 )
2523                   hom(nzb,2,77,:)       = 0.0_wp    ! because zu(nzb) is negative
2524                   data_output_pr(i)     = data_output_pr(i)(2:)
2525                ENDIF
2526             ENDIF
2527
2528          CASE ( 'w"sa"' )
2529             IF ( .NOT. ocean ) THEN
2530                message_string = 'data_output_pr = ' //                        &
2531                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2532                                 'lemented for ocean = .FALSE.'
2533                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2534             ELSE
2535                dopr_index(i) = 65
2536                dopr_unit(i)  = 'psu m/s'
2537                hom(:,2,65,:) = SPREAD( zw, 2, statistic_regions+1 )
2538             ENDIF
2539
2540          CASE ( 'w*sa*' )
2541             IF ( .NOT. ocean ) THEN
2542                message_string = 'data_output_pr = ' //                        &
2543                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2544                                 'lemented for ocean = .FALSE.'
2545                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2546             ELSE
2547                dopr_index(i) = 66
2548                dopr_unit(i)  = 'psu m/s'
2549                hom(:,2,66,:) = SPREAD( zw, 2, statistic_regions+1 )
2550             ENDIF
2551
2552          CASE ( 'wsa' )
2553             IF ( .NOT. ocean ) THEN
2554                message_string = 'data_output_pr = ' //                        &
2555                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2556                                 'lemented for ocean = .FALSE.'
2557                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2558             ELSE
2559                dopr_index(i) = 67
2560                dopr_unit(i)  = 'psu m/s'
2561                hom(:,2,67,:) = SPREAD( zw, 2, statistic_regions+1 )
2562             ENDIF
2563
2564          CASE ( 'w*p*' )
2565             dopr_index(i) = 68
2566             dopr_unit(i)  = 'm3/s3'
2567             hom(:,2,68,:) = SPREAD( zu, 2, statistic_regions+1 )
2568
2569          CASE ( 'w"e' )
2570             dopr_index(i) = 69
2571             dopr_unit(i)  = 'm3/s3'
2572             hom(:,2,69,:) = SPREAD( zu, 2, statistic_regions+1 )
2573
2574          CASE ( 'q*2' )
2575             IF ( .NOT. humidity )  THEN
2576                message_string = 'data_output_pr = ' //                        &
2577                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2578                                 'lemented for humidity = .FALSE.'
2579                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2580             ELSE
2581                dopr_index(i) = 70
2582                dopr_unit(i)  = 'kg2/kg2'
2583                hom(:,2,70,:) = SPREAD( zu, 2, statistic_regions+1 )
2584             ENDIF
2585
2586          CASE ( 'prho' )
2587             IF ( .NOT. ocean ) THEN
2588                message_string = 'data_output_pr = ' //                        &
2589                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2590                                 'lemented for ocean = .FALSE.'
2591                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2592             ELSE
2593                dopr_index(i) = 71
2594                dopr_unit(i)  = 'kg/m3'
2595                hom(:,2,71,:) = SPREAD( zu, 2, statistic_regions+1 )
2596             ENDIF
2597
2598          CASE ( 'hyp' )
2599             dopr_index(i) = 72
2600             dopr_unit(i)  = 'dbar'
2601             hom(:,2,72,:) = SPREAD( zu, 2, statistic_regions+1 )
2602
2603          CASE ( 'nr' )
2604             IF ( .NOT. cloud_physics )  THEN
2605                message_string = 'data_output_pr = ' //                        &
2606                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2607                                 'lemented for cloud_physics = .FALSE.'
2608                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
2609             ELSEIF ( icloud_scheme /= 0 )  THEN
2610                message_string = 'data_output_pr = ' //                        &
2611                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2612                                 'lemented for cloud_scheme /= seifert_beheng'
2613                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
2614             ELSEIF ( .NOT. precipitation )  THEN
2615                message_string = 'data_output_pr = ' //                        &
2616                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2617                                 'lemented for precipitation = .FALSE.'
2618                CALL message( 'check_parameters', 'PA0361', 1, 2, 0, 6, 0 )
2619             ELSE
2620                dopr_index(i) = 73
2621                dopr_unit(i)  = '1/m3'
2622                hom(:,2,73,:)  = SPREAD( zu, 2, statistic_regions+1 )
2623             ENDIF
2624
2625          CASE ( 'qr' )
2626             IF ( .NOT. cloud_physics )  THEN
2627                message_string = 'data_output_pr = ' //                        &
2628                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2629                                 'lemented for cloud_physics = .FALSE.'
2630                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
2631             ELSEIF ( icloud_scheme /= 0 )  THEN
2632                message_string = 'data_output_pr = ' //                        &
2633                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2634                                 'lemented for cloud_scheme /= seifert_beheng'
2635                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
2636             ELSEIF ( .NOT. precipitation )  THEN
2637                message_string = 'data_output_pr = ' //                        &
2638                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2639                                 'lemented for precipitation = .FALSE.'
2640                CALL message( 'check_parameters', 'PA0361', 1, 2, 0, 6, 0 )
2641             ELSE
2642                dopr_index(i) = 74
2643                dopr_unit(i)  = 'kg/kg'
2644                hom(:,2,74,:)  = SPREAD( zu, 2, statistic_regions+1 )
2645             ENDIF
2646
2647          CASE ( 'qc' )
2648             IF ( .NOT. cloud_physics )  THEN
2649                message_string = 'data_output_pr = ' //                        &
2650                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2651                                 'lemented for cloud_physics = .FALSE.'
2652                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
2653             ELSEIF ( icloud_scheme /= 0 )  THEN
2654                message_string = 'data_output_pr = ' //                        &
2655                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2656                                 'lemented for cloud_scheme /= seifert_beheng'
2657                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
2658             ELSE
2659                dopr_index(i) = 75
2660                dopr_unit(i)  = 'kg/kg'
2661                hom(:,2,75,:)  = SPREAD( zu, 2, statistic_regions+1 )
2662             ENDIF
2663
2664          CASE ( 'prr' )
2665             IF ( .NOT. cloud_physics )  THEN
2666                message_string = 'data_output_pr = ' //                        &
2667                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2668                                 'lemented for cloud_physics = .FALSE.'
2669                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
2670             ELSEIF ( icloud_scheme /= 0 )  THEN
2671                message_string = 'data_output_pr = ' //                        &
2672                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2673                                 'lemented for cloud_scheme /= seifert_beheng'
2674                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
2675             ELSEIF ( .NOT. precipitation )  THEN
2676                message_string = 'data_output_pr = ' //                        &
2677                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2678                                 'lemented for precipitation = .FALSE.'
2679                CALL message( 'check_parameters', 'PA0361', 1, 2, 0, 6, 0 )
2680
2681             ELSE
2682                dopr_index(i) = 76
2683                dopr_unit(i)  = 'kg/kg m/s'
2684                hom(:,2,76,:)  = SPREAD( zu, 2, statistic_regions+1 )
2685             ENDIF
2686
2687          CASE ( 'ug' )
2688             dopr_index(i) = 78
2689             dopr_unit(i)  = 'm/s'
2690             hom(:,2,78,:) = SPREAD( zu, 2, statistic_regions+1 )
2691
2692          CASE ( 'vg' )
2693             dopr_index(i) = 79
2694             dopr_unit(i)  = 'm/s'
2695             hom(:,2,79,:) = SPREAD( zu, 2, statistic_regions+1 )
2696
2697          CASE ( 'w_subs' )
2698             IF ( .NOT. large_scale_subsidence )  THEN
2699                message_string = 'data_output_pr = ' //                        &
2700                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2701                                 'lemented for large_scale_subsidence = .FALSE.'
2702                CALL message( 'check_parameters', 'PA0382', 1, 2, 0, 6, 0 )
2703             ELSE
2704                dopr_index(i) = 80
2705                dopr_unit(i)  = 'm/s'
2706                hom(:,2,80,:) = SPREAD( zu, 2, statistic_regions+1 )
2707             ENDIF
2708
2709          CASE ( 'td_lsa_lpt' )
2710             IF ( .NOT. large_scale_forcing )  THEN
2711                message_string = 'data_output_pr = ' //                        &
2712                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2713                                 'lemented for large_scale_forcing = .FALSE.'
2714                CALL message( 'check_parameters', 'PA0393', 1, 2, 0, 6, 0 )
2715             ELSE
2716                dopr_index(i) = 81
2717                dopr_unit(i)  = 'K/s'
2718                hom(:,2,81,:) = SPREAD( zu, 2, statistic_regions+1 )
2719             ENDIF
2720
2721          CASE ( 'td_lsa_q' )
2722             IF ( .NOT. large_scale_forcing )  THEN
2723                message_string = 'data_output_pr = ' //                        &
2724                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2725                                 'lemented for large_scale_forcing = .FALSE.'
2726                CALL message( 'check_parameters', 'PA0393', 1, 2, 0, 6, 0 )
2727             ELSE
2728                dopr_index(i) = 82
2729                dopr_unit(i)  = 'kg/kgs'
2730                hom(:,2,82,:) = SPREAD( zu, 2, statistic_regions+1 )
2731             ENDIF
2732
2733          CASE ( 'td_sub_lpt' )
2734             IF ( .NOT. large_scale_forcing )  THEN
2735                message_string = 'data_output_pr = ' //                        &
2736                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2737                                 'lemented for large_scale_forcing = .FALSE.'
2738                CALL message( 'check_parameters', 'PA0393', 1, 2, 0, 6, 0 )
2739             ELSE
2740                dopr_index(i) = 83
2741                dopr_unit(i)  = 'K/s'
2742                hom(:,2,83,:) = SPREAD( zu, 2, statistic_regions+1 )
2743             ENDIF
2744
2745          CASE ( 'td_sub_q' )
2746             IF ( .NOT. large_scale_forcing )  THEN
2747                message_string = 'data_output_pr = ' //                        &
2748                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2749                                 'lemented for large_scale_forcing = .FALSE.'
2750                CALL message( 'check_parameters', 'PA0393', 1, 2, 0, 6, 0 )
2751             ELSE
2752                dopr_index(i) = 84
2753                dopr_unit(i)  = 'kg/kgs'
2754                hom(:,2,84,:) = SPREAD( zu, 2, statistic_regions+1 )
2755             ENDIF
2756
2757          CASE ( 'td_nud_lpt' )
2758             IF ( .NOT. nudging )  THEN
2759                message_string = 'data_output_pr = ' //                        &
2760                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2761                                 'lemented for nudging = .FALSE.'
2762                CALL message( 'check_parameters', 'PA0394', 1, 2, 0, 6, 0 )
2763             ELSE
2764                dopr_index(i) = 85
2765                dopr_unit(i)  = 'K/s'
2766                hom(:,2,85,:) = SPREAD( zu, 2, statistic_regions+1 )
2767             ENDIF
2768
2769          CASE ( 'td_nud_q' )
2770             IF ( .NOT. nudging )  THEN
2771                message_string = 'data_output_pr = ' //                        &
2772                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2773                                 'lemented for nudging = .FALSE.'
2774                CALL message( 'check_parameters', 'PA0394', 1, 2, 0, 6, 0 )
2775             ELSE
2776                dopr_index(i) = 86
2777                dopr_unit(i)  = 'kg/kgs'
2778                hom(:,2,86,:) = SPREAD( zu, 2, statistic_regions+1 )
2779             ENDIF
2780
2781          CASE ( 'td_nud_u' )
2782             IF ( .NOT. nudging )  THEN
2783                message_string = 'data_output_pr = ' //                        &
2784                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2785                                 'lemented for nudging = .FALSE.'
2786                CALL message( 'check_parameters', 'PA0394', 1, 2, 0, 6, 0 )
2787             ELSE
2788                dopr_index(i) = 87
2789                dopr_unit(i)  = 'm/s2'
2790                hom(:,2,87,:) = SPREAD( zu, 2, statistic_regions+1 )
2791             ENDIF
2792
2793          CASE ( 'td_nud_v' )
2794             IF ( .NOT. nudging )  THEN
2795                message_string = 'data_output_pr = ' //                        &
2796                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2797                                 'lemented for nudging = .FALSE.'
2798                CALL message( 'check_parameters', 'PA0394', 1, 2, 0, 6, 0 )
2799             ELSE
2800                dopr_index(i) = 88
2801                dopr_unit(i)  = 'm/s2'
2802                hom(:,2,88,:) = SPREAD( zu, 2, statistic_regions+1 )
2803             ENDIF
2804
2805
2806          CASE DEFAULT
2807
2808             CALL user_check_data_output_pr( data_output_pr(i), i, unit )
2809
2810             IF ( unit == 'illegal' )  THEN
2811                IF ( data_output_pr_user(1) /= ' ' )  THEN
2812                   message_string = 'illegal value for data_output_pr or ' //  &
2813                                    'data_output_pr_user = "' //               &
2814                                    TRIM( data_output_pr(i) ) // '"'
2815                   CALL message( 'check_parameters', 'PA0097', 1, 2, 0, 6, 0 )
2816                ELSE
2817                   message_string = 'illegal value for data_output_pr = "' //  &
2818                                    TRIM( data_output_pr(i) ) // '"'
2819                   CALL message( 'check_parameters', 'PA0098', 1, 2, 0, 6, 0 )
2820                ENDIF
2821             ENDIF
2822
2823       END SELECT
2824
2825    ENDDO
2826
2827
2828!
2829!-- Append user-defined data output variables to the standard data output
2830    IF ( data_output_user(1) /= ' ' )  THEN
2831       i = 1
2832       DO  WHILE ( data_output(i) /= ' '  .AND.  i <= 100 )
2833          i = i + 1
2834       ENDDO
2835       j = 1
2836       DO  WHILE ( data_output_user(j) /= ' '  .AND.  j <= 100 )
2837          IF ( i > 100 )  THEN
2838             message_string = 'number of output quantitities given by data' // &
2839                '_output and data_output_user exceeds the limit of 100'
2840             CALL message( 'check_parameters', 'PA0102', 1, 2, 0, 6, 0 )
2841          ENDIF
2842          data_output(i) = data_output_user(j)
2843          i = i + 1
2844          j = j + 1
2845       ENDDO
2846    ENDIF
2847
2848!
2849!-- Check and set steering parameters for 2d/3d data output and averaging
2850    i   = 1
2851    DO  WHILE ( data_output(i) /= ' '  .AND.  i <= 100 )
2852!
2853!--    Check for data averaging
2854       ilen = LEN_TRIM( data_output(i) )
2855       j = 0                                                 ! no data averaging
2856       IF ( ilen > 3 )  THEN
2857          IF ( data_output(i)(ilen-2:ilen) == '_av' )  THEN
2858             j = 1                                           ! data averaging
2859             data_output(i) = data_output(i)(1:ilen-3)
2860          ENDIF
2861       ENDIF
2862!
2863!--    Check for cross section or volume data
2864       ilen = LEN_TRIM( data_output(i) )
2865       k = 0                                                   ! 3d data
2866       var = data_output(i)(1:ilen)
2867       IF ( ilen > 3 )  THEN
2868          IF ( data_output(i)(ilen-2:ilen) == '_xy'  .OR.                      &
2869               data_output(i)(ilen-2:ilen) == '_xz'  .OR.                      &
2870               data_output(i)(ilen-2:ilen) == '_yz' )  THEN
2871             k = 1                                             ! 2d data
2872             var = data_output(i)(1:ilen-3)
2873          ENDIF
2874       ENDIF
2875!
2876!--    Check for allowed value and set units
2877       SELECT CASE ( TRIM( var ) )
2878
2879          CASE ( 'e' )
2880             IF ( constant_diffusion )  THEN
2881                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
2882                                 'res constant_diffusion = .FALSE.'
2883                CALL message( 'check_parameters', 'PA0103', 1, 2, 0, 6, 0 )
2884             ENDIF
2885             unit = 'm2/s2'
2886
2887          CASE ( 'lpt' )
2888             IF ( .NOT. cloud_physics )  THEN
2889                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
2890                         'res cloud_physics = .TRUE.'
2891                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
2892             ENDIF
2893             unit = 'K'
2894
2895          CASE ( 'nr' )
2896             IF ( .NOT. cloud_physics )  THEN
2897                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
2898                         'res cloud_physics = .TRUE.'
2899                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
2900             ELSEIF ( icloud_scheme /= 0 )  THEN
2901                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
2902                         'res cloud_scheme = seifert_beheng'
2903                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
2904             ENDIF
2905             unit = '1/m3'
2906
2907          CASE ( 'pc', 'pr' )
2908             IF ( .NOT. particle_advection )  THEN
2909                message_string = 'output of "' // TRIM( var ) // '" requir' // &
2910                   'es a "particles_par"-NAMELIST in the parameter file (PARIN)'
2911                CALL message( 'check_parameters', 'PA0104', 1, 2, 0, 6, 0 )
2912             ENDIF
2913             IF ( TRIM( var ) == 'pc' )  unit = 'number'
2914             IF ( TRIM( var ) == 'pr' )  unit = 'm'
2915
2916          CASE ( 'prr' )
2917             IF ( .NOT. cloud_physics )  THEN
2918                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
2919                         'res cloud_physics = .TRUE.'
2920                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
2921             ELSEIF ( icloud_scheme /= 0 )  THEN
2922                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
2923                         'res cloud_scheme = seifert_beheng'
2924                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
2925             ELSEIF ( .NOT. precipitation )  THEN
2926                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
2927                                 'res precipitation = .TRUE.'
2928                CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
2929             ENDIF
2930             unit = 'kg/kg m/s'
2931
2932          CASE ( 'q', 'vpt' )
2933             IF ( .NOT. humidity )  THEN
2934                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
2935                                 'res humidity = .TRUE.'
2936                CALL message( 'check_parameters', 'PA0105', 1, 2, 0, 6, 0 )
2937             ENDIF
2938             IF ( TRIM( var ) == 'q'   )  unit = 'kg/kg'
2939             IF ( TRIM( var ) == 'vpt' )  unit = 'K'
2940
2941          CASE ( 'qc' )
2942             IF ( .NOT. cloud_physics )  THEN
2943                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
2944                         'res cloud_physics = .TRUE.'
2945                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
2946             ELSEIF ( icloud_scheme /= 0 ) THEN
2947                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
2948                         'res cloud_scheme = seifert_beheng'
2949                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
2950             ENDIF
2951             unit = 'kg/kg'
2952
2953          CASE ( 'ql' )
2954             IF ( .NOT. ( cloud_physics  .OR.  cloud_droplets ) )  THEN
2955                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
2956                         'res cloud_physics = .TRUE. or cloud_droplets = .TRUE.'
2957                CALL message( 'check_parameters', 'PA0106', 1, 2, 0, 6, 0 )
2958             ENDIF
2959             unit = 'kg/kg'
2960
2961          CASE ( 'ql_c', 'ql_v', 'ql_vp' )
2962             IF ( .NOT. cloud_droplets )  THEN
2963                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
2964                                 'res cloud_droplets = .TRUE.'
2965                CALL message( 'check_parameters', 'PA0107', 1, 2, 0, 6, 0 )
2966             ENDIF
2967             IF ( TRIM( var ) == 'ql_c'  )  unit = 'kg/kg'
2968             IF ( TRIM( var ) == 'ql_v'  )  unit = 'm3'
2969             IF ( TRIM( var ) == 'ql_vp' )  unit = 'none'
2970
2971          CASE ( 'qr' )
2972             IF ( .NOT. cloud_physics )  THEN
2973                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
2974                         'res cloud_physics = .TRUE.'
2975                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
2976             ELSEIF ( icloud_scheme /= 0 ) THEN
2977                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
2978                         'res cloud_scheme = seifert_beheng'
2979                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
2980             ELSEIF ( .NOT. precipitation )  THEN
2981                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
2982                                 'res precipitation = .TRUE.'
2983                CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
2984             ENDIF
2985             unit = 'kg/kg'
2986
2987          CASE ( 'qv' )
2988             IF ( .NOT. cloud_physics )  THEN
2989                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
2990                                 'res cloud_physics = .TRUE.'
2991                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
2992             ENDIF
2993             unit = 'kg/kg'
2994
2995          CASE ( 'rho' )
2996             IF ( .NOT. ocean )  THEN
2997                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
2998                                 'res ocean = .TRUE.'
2999                CALL message( 'check_parameters', 'PA0109', 1, 2, 0, 6, 0 )
3000             ENDIF
3001             unit = 'kg/m3'
3002
3003          CASE ( 's' )
3004             IF ( .NOT. passive_scalar )  THEN
3005                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3006                                 'res passive_scalar = .TRUE.'
3007                CALL message( 'check_parameters', 'PA0110', 1, 2, 0, 6, 0 )
3008             ENDIF
3009             unit = 'conc'
3010
3011          CASE ( 'sa' )
3012             IF ( .NOT. ocean )  THEN
3013                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3014                                 'res ocean = .TRUE.'
3015                CALL message( 'check_parameters', 'PA0109', 1, 2, 0, 6, 0 )
3016             ENDIF
3017             unit = 'psu'
3018
3019          CASE ( 'u*', 't*', 'lwp*', 'pra*', 'prr*', 'qsws*', 'shf*', 'z0*', 'z0h*' )
3020             IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
3021                message_string = 'illegal value for data_output: "' //         &
3022                                 TRIM( var ) // '" & only 2d-horizontal ' //   &
3023                                 'cross sections are allowed for this value'
3024                CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
3025             ENDIF
3026             IF ( TRIM( var ) == 'lwp*'  .AND.  .NOT. cloud_physics )  THEN
3027                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3028                                 'res cloud_physics = .TRUE.'
3029                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3030             ENDIF
3031             IF ( TRIM( var ) == 'pra*'  .AND.  .NOT. precipitation )  THEN
3032                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3033                                 'res precipitation = .TRUE.'
3034                CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
3035             ENDIF
3036             IF ( TRIM( var ) == 'pra*'  .AND.  j == 1 )  THEN
3037                message_string = 'temporal averaging of precipitation ' //     &
3038                          'amount "' // TRIM( var ) // '" is not possible'
3039                CALL message( 'check_parameters', 'PA0113', 1, 2, 0, 6, 0 )
3040             ENDIF
3041             IF ( TRIM( var ) == 'prr*'  .AND.  .NOT. precipitation )  THEN
3042                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3043                                 'res precipitation = .TRUE.'
3044                CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
3045             ENDIF
3046             IF ( TRIM( var ) == 'qsws*'  .AND.  .NOT. humidity )  THEN
3047                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3048                                 'res humidity = .TRUE.'
3049                CALL message( 'check_parameters', 'PA0322', 1, 2, 0, 6, 0 )
3050             ENDIF
3051
3052             IF ( TRIM( var ) == 'lwp*'   )  unit = 'kg/kg*m'
3053             IF ( TRIM( var ) == 'pra*'   )  unit = 'mm'
3054             IF ( TRIM( var ) == 'prr*'   )  unit = 'mm/s'
3055             IF ( TRIM( var ) == 'qsws*'  )  unit = 'kgm/kgs'
3056             IF ( TRIM( var ) == 'shf*'   )  unit = 'K*m/s'
3057             IF ( TRIM( var ) == 't*'     )  unit = 'K'
3058             IF ( TRIM( var ) == 'u*'     )  unit = 'm/s'
3059             IF ( TRIM( var ) == 'z0*'    )  unit = 'm'
3060             IF ( TRIM( var ) == 'z0h*'   )  unit = 'm'
3061
3062
3063          CASE ( 'p', 'pt', 'u', 'v', 'w' )
3064             IF ( TRIM( var ) == 'p'  )  unit = 'Pa'
3065             IF ( TRIM( var ) == 'pt' )  unit = 'K'
3066             IF ( TRIM( var ) == 'u'  )  unit = 'm/s'
3067             IF ( TRIM( var ) == 'v'  )  unit = 'm/s'
3068             IF ( TRIM( var ) == 'w'  )  unit = 'm/s'
3069             CONTINUE
3070
3071          CASE DEFAULT
3072             CALL user_check_data_output( var, unit )
3073
3074             IF ( unit == 'illegal' )  THEN
3075                IF ( data_output_user(1) /= ' ' )  THEN
3076                   message_string = 'illegal value for data_output or ' //     &
3077                         'data_output_user = "' // TRIM( data_output(i) ) // '"'
3078                   CALL message( 'check_parameters', 'PA0114', 1, 2, 0, 6, 0 )
3079                ELSE
3080                   message_string = 'illegal value for data_output =' //       &
3081                                    TRIM( data_output(i) ) // '"'
3082                   CALL message( 'check_parameters', 'PA0115', 1, 2, 0, 6, 0 )
3083                ENDIF
3084             ENDIF
3085
3086       END SELECT
3087!
3088!--    Set the internal steering parameters appropriately
3089       IF ( k == 0 )  THEN
3090          do3d_no(j)              = do3d_no(j) + 1
3091          do3d(j,do3d_no(j))      = data_output(i)
3092          do3d_unit(j,do3d_no(j)) = unit
3093       ELSE
3094          do2d_no(j)              = do2d_no(j) + 1
3095          do2d(j,do2d_no(j))      = data_output(i)
3096          do2d_unit(j,do2d_no(j)) = unit
3097          IF ( data_output(i)(ilen-2:ilen) == '_xy' )  THEN
3098             data_output_xy(j) = .TRUE.
3099          ENDIF
3100          IF ( data_output(i)(ilen-2:ilen) == '_xz' )  THEN
3101             data_output_xz(j) = .TRUE.
3102          ENDIF
3103          IF ( data_output(i)(ilen-2:ilen) == '_yz' )  THEN
3104             data_output_yz(j) = .TRUE.
3105          ENDIF
3106       ENDIF
3107
3108       IF ( j == 1 )  THEN
3109!
3110!--       Check, if variable is already subject to averaging
3111          found = .FALSE.
3112          DO  k = 1, doav_n
3113             IF ( TRIM( doav(k) ) == TRIM( var ) )  found = .TRUE.
3114          ENDDO
3115
3116          IF ( .NOT. found )  THEN
3117             doav_n = doav_n + 1
3118             doav(doav_n) = var
3119          ENDIF
3120       ENDIF
3121
3122       i = i + 1
3123    ENDDO
3124
3125!
3126!-- Averaged 2d or 3d output requires that an averaging interval has been set
3127    IF ( doav_n > 0  .AND.  averaging_interval == 0.0_wp )  THEN
3128       WRITE( message_string, * )  'output of averaged quantity "',            &
3129                                   TRIM( doav(1) ), '_av" requires to set a ', &
3130                                   'non-zero & averaging interval'
3131       CALL message( 'check_parameters', 'PA0323', 1, 2, 0, 6, 0 )
3132    ENDIF
3133
3134!
3135!-- Check sectional planes and store them in one shared array
3136    IF ( ANY( section_xy > nz + 1 ) )  THEN
3137       WRITE( message_string, * )  'section_xy must be <= nz + 1 = ', nz + 1
3138       CALL message( 'check_parameters', 'PA0319', 1, 2, 0, 6, 0 )
3139    ENDIF
3140    IF ( ANY( section_xz > ny + 1 ) )  THEN
3141       WRITE( message_string, * )  'section_xz must be <= ny + 1 = ', ny + 1
3142       CALL message( 'check_parameters', 'PA0320', 1, 2, 0, 6, 0 )
3143    ENDIF
3144    IF ( ANY( section_yz > nx + 1 ) )  THEN
3145       WRITE( message_string, * )  'section_yz must be <= nx + 1 = ', nx + 1
3146       CALL message( 'check_parameters', 'PA0321', 1, 2, 0, 6, 0 )
3147    ENDIF
3148    section(:,1) = section_xy
3149    section(:,2) = section_xz
3150    section(:,3) = section_yz
3151
3152!
3153!-- Upper plot limit for 2D vertical sections
3154    IF ( z_max_do2d == -1.0_wp )  z_max_do2d = zu(nzt)
3155    IF ( z_max_do2d < zu(nzb+1)  .OR.  z_max_do2d > zu(nzt) )  THEN
3156       WRITE( message_string, * )  'z_max_do2d = ', z_max_do2d,                &
3157                    ' must be >= ', zu(nzb+1), '(zu(nzb+1)) and <= ', zu(nzt), &
3158                    ' (zu(nzt))'
3159       CALL message( 'check_parameters', 'PA0116', 1, 2, 0, 6, 0 )
3160    ENDIF
3161
3162!
3163!-- Upper plot limit for 3D arrays
3164    IF ( nz_do3d == -9999 )  nz_do3d = nzt + 1
3165
3166!
3167!-- Set output format string (used in header)
3168    SELECT CASE ( netcdf_data_format )
3169       CASE ( 1 )
3170          output_format_netcdf = 'netCDF classic'
3171       CASE ( 2 )
3172          output_format_netcdf = 'netCDF 64bit offset'
3173       CASE ( 3 )
3174          output_format_netcdf = 'netCDF4/HDF5'
3175       CASE ( 4 )
3176          output_format_netcdf = 'netCDF4/HDF5 classic'
3177       CASE ( 5 )
3178          output_format_netcdf = 'parallel netCDF4/HDF5'
3179       CASE ( 6 )
3180          output_format_netcdf = 'parallel netCDF4/HDF5 classic'
3181
3182    END SELECT
3183
3184#if defined( __spectra )
3185!
3186!-- Check the number of spectra level to be output
3187    i = 1
3188    DO WHILE ( comp_spectra_level(i) /= 999999  .AND.  i <= 100 )
3189       i = i + 1
3190    ENDDO
3191    i = i - 1
3192    IF ( i == 0 )  THEN
3193       WRITE( message_string, * )  'no spectra levels given'
3194       CALL message( 'check_parameters', 'PA0019', 1, 2, 0, 6, 0 )
3195    ENDIF
3196#endif
3197
3198!
3199!-- Check mask conditions
3200    DO mid = 1, max_masks
3201       IF ( data_output_masks(mid,1) /= ' ' .OR.                               &
3202            data_output_masks_user(mid,1) /= ' ' ) THEN
3203          masks = masks + 1
3204       ENDIF
3205    ENDDO
3206   
3207    IF ( masks < 0 .OR. masks > max_masks )  THEN
3208       WRITE( message_string, * )  'illegal value: masks must be >= 0 and ',   &
3209            '<= ', max_masks, ' (=max_masks)'
3210       CALL message( 'check_parameters', 'PA0325', 1, 2, 0, 6, 0 )
3211    ENDIF
3212    IF ( masks > 0 )  THEN
3213       mask_scale(1) = mask_scale_x
3214       mask_scale(2) = mask_scale_y
3215       mask_scale(3) = mask_scale_z
3216       IF ( ANY( mask_scale <= 0.0_wp ) )  THEN
3217          WRITE( message_string, * )                                           &
3218               'illegal value: mask_scale_x, mask_scale_y and mask_scale_z',   &
3219               'must be > 0.0'
3220          CALL message( 'check_parameters', 'PA0326', 1, 2, 0, 6, 0 )
3221       ENDIF
3222!
3223!--    Generate masks for masked data output
3224!--    Parallel netcdf output is not tested so far for masked data, hence
3225!--    netcdf_data_format is switched back to non-paralell output.
3226       netcdf_data_format_save = netcdf_data_format
3227       IF ( netcdf_data_format > 4 )  THEN
3228          IF ( netcdf_data_format == 5 ) netcdf_data_format = 3
3229          IF ( netcdf_data_format == 6 ) netcdf_data_format = 4
3230          message_string = 'netCDF file formats '//                            &
3231                           '5 (parallel netCDF 4) and ' //                     &
3232                           '6 (parallel netCDF 4 Classic model) '//            &
3233                           '&are currently not supported (not yet tested) ' // &
3234                           'for masked data.&Using respective non-parallel' // & 
3235                           ' output for masked data.'
3236          CALL message( 'check_parameters', 'PA0383', 0, 0, 0, 6, 0 )
3237       ENDIF
3238       CALL init_masks
3239       netcdf_data_format = netcdf_data_format_save
3240    ENDIF
3241
3242!
3243!-- Check the NetCDF data format
3244#if ! defined ( __check )
3245    IF ( netcdf_data_format > 2 )  THEN
3246#if defined( __netcdf4 )
3247       CONTINUE
3248#else
3249       message_string = 'netCDF: netCDF4 format requested but no ' //          &
3250                        'cpp-directive __netcdf4 given & switch '  //          &
3251                        'back to 64-bit offset format'
3252       CALL message( 'check_parameters', 'PA0171', 0, 1, 0, 6, 0 )
3253       netcdf_data_format = 2
3254#endif
3255    ENDIF
3256    IF ( netcdf_data_format > 4 )  THEN
3257#if defined( __netcdf4 ) && defined( __netcdf4_parallel )
3258       CONTINUE
3259#else
3260       message_string = 'netCDF: netCDF4 parallel output requested but no ' // &
3261                        'cpp-directive __netcdf4_parallel given & switch '  // &
3262                        'back to netCDF4 non-parallel output'
3263       CALL message( 'check_parameters', 'PA0099', 0, 1, 0, 6, 0 )
3264       netcdf_data_format = netcdf_data_format - 2
3265#endif
3266    ENDIF
3267#endif
3268
3269!
3270!-- Calculate fixed number of output time levels for parallel netcdf output.
3271!-- The time dimension has to be defined as limited for parallel output,
3272!-- because otherwise the I/O performance drops significantly.
3273    IF ( netcdf_data_format > 4 )  THEN
3274
3275       ntdim_3d(0)    = INT( ( end_time - skip_time_do3d ) / dt_do3d )
3276       IF ( do3d_at_begin ) ntdim_3d(0) = ntdim_3d(0) + 1
3277       ntdim_3d(1)    = INT( ( end_time - skip_time_data_output_av )           &
3278                             / dt_data_output_av )
3279       ntdim_2d_xy(0) = INT( ( end_time - skip_time_do2d_xy ) / dt_do2d_xy )
3280       ntdim_2d_xz(0) = INT( ( end_time - skip_time_do2d_xz ) / dt_do2d_xz )
3281       ntdim_2d_yz(0) = INT( ( end_time - skip_time_do2d_yz ) / dt_do2d_yz )
3282       IF ( do2d_at_begin )  THEN
3283          ntdim_2d_xy(0) = ntdim_2d_xy(0) + 1
3284          ntdim_2d_xz(0) = ntdim_2d_xz(0) + 1
3285          ntdim_2d_yz(0) = ntdim_2d_yz(0) + 1
3286       ENDIF
3287       ntdim_2d_xy(1) = ntdim_3d(1)
3288       ntdim_2d_xz(1) = ntdim_3d(1)
3289       ntdim_2d_yz(1) = ntdim_3d(1)
3290
3291    ENDIF
3292
3293#if ! defined( __check )
3294!
3295!-- Check netcdf precison
3296    ldum = .FALSE.
3297    CALL define_netcdf_header( 'ch', ldum, 0 )
3298#endif
3299!
3300!-- Check, whether a constant diffusion coefficient shall be used
3301    IF ( km_constant /= -1.0_wp )  THEN
3302       IF ( km_constant < 0.0_wp )  THEN
3303          WRITE( message_string, * )  'km_constant = ', km_constant, ' < 0.0'
3304          CALL message( 'check_parameters', 'PA0121', 1, 2, 0, 6, 0 )
3305       ELSE
3306          IF ( prandtl_number < 0.0_wp )  THEN
3307             WRITE( message_string, * )  'prandtl_number = ', prandtl_number,  &
3308                                         ' < 0.0'
3309             CALL message( 'check_parameters', 'PA0122', 1, 2, 0, 6, 0 )
3310          ENDIF
3311          constant_diffusion = .TRUE.
3312
3313          IF ( prandtl_layer )  THEN
3314             message_string = 'prandtl_layer is not allowed with fixed ' //    &
3315                              'value of km'
3316             CALL message( 'check_parameters', 'PA0123', 1, 2, 0, 6, 0 )
3317          ENDIF
3318       ENDIF
3319    ENDIF
3320
3321!
3322!-- In case of non-cyclic lateral boundaries and a damping layer for the
3323!-- potential temperature, check the width of the damping layer
3324    IF ( bc_lr /= 'cyclic' ) THEN
3325       IF ( pt_damping_width < 0.0_wp  .OR.                                    &
3326            pt_damping_width > REAL( nx * dx ) )  THEN
3327          message_string = 'pt_damping_width out of range'
3328          CALL message( 'check_parameters', 'PA0124', 1, 2, 0, 6, 0 )
3329       ENDIF
3330    ENDIF
3331
3332    IF ( bc_ns /= 'cyclic' )  THEN
3333       IF ( pt_damping_width < 0.0_wp  .OR.                                    &
3334            pt_damping_width > REAL( ny * dy ) )  THEN
3335          message_string = 'pt_damping_width out of range'
3336          CALL message( 'check_parameters', 'PA0124', 1, 2, 0, 6, 0 )
3337       ENDIF
3338    ENDIF
3339
3340!
3341!-- Check value range for rif
3342    IF ( rif_min >= rif_max )  THEN
3343       WRITE( message_string, * )  'rif_min = ', rif_min, ' must be less ',    &
3344                                   'than rif_max = ', rif_max
3345       CALL message( 'check_parameters', 'PA0125', 1, 2, 0, 6, 0 )
3346    ENDIF
3347
3348!
3349!-- Check random generator
3350    IF ( (random_generator /= 'system-specific'     .AND.                      &
3351          random_generator /= 'random-parallel'   ) .AND.                      &
3352          random_generator /= 'numerical-recipes' )  THEN
3353       message_string = 'unknown random generator: random_generator = "' //    &
3354                        TRIM( random_generator ) // '"'
3355       CALL message( 'check_parameters', 'PA0135', 1, 2, 0, 6, 0 )
3356    ENDIF
3357
3358!
3359!-- Determine upper and lower hight level indices for random perturbations
3360    IF ( disturbance_level_b == -9999999.9_wp )  THEN
3361       IF ( ocean ) THEN
3362          disturbance_level_b     = zu((nzt*2)/3)
3363          disturbance_level_ind_b = ( nzt * 2 ) / 3
3364       ELSE
3365          disturbance_level_b     = zu(nzb+3)
3366          disturbance_level_ind_b = nzb + 3
3367       ENDIF
3368    ELSEIF ( disturbance_level_b < zu(3) )  THEN
3369       WRITE( message_string, * )  'disturbance_level_b = ',                   &
3370                           disturbance_level_b, ' must be >= ', zu(3), '(zu(3))'
3371       CALL message( 'check_parameters', 'PA0126', 1, 2, 0, 6, 0 )
3372    ELSEIF ( disturbance_level_b > zu(nzt-2) )  THEN
3373       WRITE( message_string, * )  'disturbance_level_b = ',                   &
3374                   disturbance_level_b, ' must be <= ', zu(nzt-2), '(zu(nzt-2))'
3375       CALL message( 'check_parameters', 'PA0127', 1, 2, 0, 6, 0 )
3376    ELSE
3377       DO  k = 3, nzt-2
3378          IF ( disturbance_level_b <= zu(k) )  THEN
3379             disturbance_level_ind_b = k
3380             EXIT
3381          ENDIF
3382       ENDDO
3383    ENDIF
3384
3385    IF ( disturbance_level_t == -9999999.9_wp )  THEN
3386       IF ( ocean )  THEN
3387          disturbance_level_t     = zu(nzt-3)
3388          disturbance_level_ind_t = nzt - 3
3389       ELSE
3390          disturbance_level_t     = zu(nzt/3)
3391          disturbance_level_ind_t = nzt / 3
3392       ENDIF
3393    ELSEIF ( disturbance_level_t > zu(nzt-2) )  THEN
3394       WRITE( message_string, * )  'disturbance_level_t = ',                   &
3395                   disturbance_level_t, ' must be <= ', zu(nzt-2), '(zu(nzt-2))'
3396       CALL message( 'check_parameters', 'PA0128', 1, 2, 0, 6, 0 )
3397    ELSEIF ( disturbance_level_t < disturbance_level_b )  THEN
3398       WRITE( message_string, * )  'disturbance_level_t = ',                   &
3399                   disturbance_level_t, ' must be >= disturbance_level_b = ',  &
3400                   disturbance_level_b
3401       CALL message( 'check_parameters', 'PA0129', 1, 2, 0, 6, 0 )
3402    ELSE
3403       DO  k = 3, nzt-2
3404          IF ( disturbance_level_t <= zu(k) )  THEN
3405             disturbance_level_ind_t = k
3406             EXIT
3407          ENDIF
3408       ENDDO
3409    ENDIF
3410
3411!
3412!-- Check again whether the levels determined this way are ok.
3413!-- Error may occur at automatic determination and too few grid points in
3414!-- z-direction.
3415    IF ( disturbance_level_ind_t < disturbance_level_ind_b )  THEN
3416       WRITE( message_string, * )  'disturbance_level_ind_t = ',               &
3417                disturbance_level_ind_t, ' must be >= disturbance_level_b = ', &
3418                disturbance_level_b
3419       CALL message( 'check_parameters', 'PA0130', 1, 2, 0, 6, 0 )
3420    ENDIF
3421
3422!
3423!-- Determine the horizontal index range for random perturbations.
3424!-- In case of non-cyclic horizontal boundaries, no perturbations are imposed
3425!-- near the inflow and the perturbation area is further limited to ...(1)
3426!-- after the initial phase of the flow.
3427   
3428    IF ( bc_lr /= 'cyclic' )  THEN
3429       IF ( inflow_disturbance_begin == -1 )  THEN
3430          inflow_disturbance_begin = MIN( 10, nx/2 )
3431       ENDIF
3432       IF ( inflow_disturbance_begin < 0  .OR.  inflow_disturbance_begin > nx )&
3433       THEN
3434          message_string = 'inflow_disturbance_begin out of range'
3435          CALL message( 'check_parameters', 'PA0131', 1, 2, 0, 6, 0 )
3436       ENDIF
3437       IF ( inflow_disturbance_end == -1 )  THEN
3438          inflow_disturbance_end = MIN( 100, 3*nx/4 )
3439       ENDIF
3440       IF ( inflow_disturbance_end < 0  .OR.  inflow_disturbance_end > nx )    &
3441       THEN
3442          message_string = 'inflow_disturbance_end out of range'
3443          CALL message( 'check_parameters', 'PA0132', 1, 2, 0, 6, 0 )
3444       ENDIF
3445    ELSEIF ( bc_ns /= 'cyclic' )  THEN
3446       IF ( inflow_disturbance_begin == -1 )  THEN
3447          inflow_disturbance_begin = MIN( 10, ny/2 )
3448       ENDIF
3449       IF ( inflow_disturbance_begin < 0  .OR.  inflow_disturbance_begin > ny )&
3450       THEN
3451          message_string = 'inflow_disturbance_begin out of range'
3452          CALL message( 'check_parameters', 'PA0131', 1, 2, 0, 6, 0 )
3453       ENDIF
3454       IF ( inflow_disturbance_end == -1 )  THEN
3455          inflow_disturbance_end = MIN( 100, 3*ny/4 )
3456       ENDIF
3457       IF ( inflow_disturbance_end < 0  .OR.  inflow_disturbance_end > ny )    &
3458       THEN
3459          message_string = 'inflow_disturbance_end out of range'
3460          CALL message( 'check_parameters', 'PA0132', 1, 2, 0, 6, 0 )
3461       ENDIF
3462    ENDIF
3463
3464    IF ( random_generator == 'random-parallel' )  THEN
3465       dist_nxl = nxl;  dist_nxr = nxr
3466       dist_nys = nys;  dist_nyn = nyn
3467       IF ( bc_lr == 'radiation/dirichlet' )  THEN
3468          dist_nxr    = MIN( nx - inflow_disturbance_begin, nxr )
3469          dist_nxl(1) = MAX( nx - inflow_disturbance_end, nxl )
3470       ELSEIF ( bc_lr == 'dirichlet/radiation' )  THEN
3471          dist_nxl    = MAX( inflow_disturbance_begin, nxl )
3472          dist_nxr(1) = MIN( inflow_disturbance_end, nxr )
3473       ENDIF
3474       IF ( bc_ns == 'dirichlet/radiation' )  THEN
3475          dist_nyn    = MIN( ny - inflow_disturbance_begin, nyn )
3476          dist_nys(1) = MAX( ny - inflow_disturbance_end, nys )
3477       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
3478          dist_nys    = MAX( inflow_disturbance_begin, nys )
3479          dist_nyn(1) = MIN( inflow_disturbance_end, nyn )
3480       ENDIF
3481    ELSE
3482       dist_nxl = 0;  dist_nxr = nx
3483       dist_nys = 0;  dist_nyn = ny
3484       IF ( bc_lr == 'radiation/dirichlet' )  THEN
3485          dist_nxr    = nx - inflow_disturbance_begin
3486          dist_nxl(1) = nx - inflow_disturbance_end
3487       ELSEIF ( bc_lr == 'dirichlet/radiation' )  THEN
3488          dist_nxl    = inflow_disturbance_begin
3489          dist_nxr(1) = inflow_disturbance_end
3490       ENDIF
3491       IF ( bc_ns == 'dirichlet/radiation' )  THEN
3492          dist_nyn    = ny - inflow_disturbance_begin
3493          dist_nys(1) = ny - inflow_disturbance_end
3494       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
3495          dist_nys    = inflow_disturbance_begin
3496          dist_nyn(1) = inflow_disturbance_end
3497       ENDIF
3498    ENDIF
3499
3500!
3501!-- A turbulent inflow requires Dirichlet conditions at the respective inflow
3502!-- boundary (so far, a turbulent inflow is realized from the left side only)
3503    IF ( turbulent_inflow  .AND.  bc_lr /= 'dirichlet/radiation' )  THEN
3504       message_string = 'turbulent_inflow = .T. requires a Dirichlet ' //      &
3505                        'condition at the inflow boundary'
3506       CALL message( 'check_parameters', 'PA0133', 1, 2, 0, 6, 0 )
3507    ENDIF
3508
3509!
3510!-- Turbulent inflow requires that 3d arrays have been cyclically filled with
3511!-- data from prerun in the first main run
3512    IF ( turbulent_inflow  .AND.  initializing_actions /= 'cyclic_fill'  .AND. &
3513         initializing_actions /= 'read_restart_data' )  THEN
3514       message_string = 'turbulent_inflow = .T. requires ' //                  &
3515                        'initializing_actions = ''cyclic_fill'' '
3516       CALL message( 'check_parameters', 'PA0055', 1, 2, 0, 6, 0 )
3517    ENDIF
3518
3519!
3520!-- In case of turbulent inflow calculate the index of the recycling plane
3521    IF ( turbulent_inflow )  THEN
3522       IF ( recycling_width == 9999999.9_wp )  THEN
3523!
3524!--       Set the default value for the width of the recycling domain
3525          recycling_width = 0.1_wp * nx * dx
3526       ELSE
3527          IF ( recycling_width < dx  .OR.  recycling_width > nx * dx )  THEN
3528             WRITE( message_string, * )  'illegal value for recycling_width:', &
3529                                         ' ', recycling_width
3530             CALL message( 'check_parameters', 'PA0134', 1, 2, 0, 6, 0 )
3531          ENDIF
3532       ENDIF
3533!
3534!--    Calculate the index
3535       recycling_plane = recycling_width / dx
3536    ENDIF
3537
3538!
3539!-- Determine damping level index for 1D model
3540    IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
3541       IF ( damp_level_1d == -1.0_wp )  THEN
3542          damp_level_1d     = zu(nzt+1)
3543          damp_level_ind_1d = nzt + 1
3544       ELSEIF ( damp_level_1d < 0.0_wp  .OR.  damp_level_1d > zu(nzt+1) )  THEN
3545          WRITE( message_string, * )  'damp_level_1d = ', damp_level_1d,       &
3546                 ' must be > 0.0 and < ', zu(nzt+1), '(zu(nzt+1))'
3547          CALL message( 'check_parameters', 'PA0136', 1, 2, 0, 6, 0 )
3548       ELSE
3549          DO  k = 1, nzt+1
3550             IF ( damp_level_1d <= zu(k) )  THEN
3551                damp_level_ind_1d = k
3552                EXIT
3553             ENDIF
3554          ENDDO
3555       ENDIF
3556    ENDIF
3557
3558!
3559!-- Check some other 1d-model parameters
3560    IF ( TRIM( mixing_length_1d ) /= 'as_in_3d_model'  .AND.                   &
3561         TRIM( mixing_length_1d ) /= 'blackadar' )  THEN
3562       message_string = 'mixing_length_1d = "' // TRIM( mixing_length_1d ) //  &
3563                        '" is unknown'
3564       CALL message( 'check_parameters', 'PA0137', 1, 2, 0, 6, 0 )
3565    ENDIF
3566    IF ( TRIM( dissipation_1d ) /= 'as_in_3d_model'  .AND.                     &
3567         TRIM( dissipation_1d ) /= 'detering' )  THEN
3568       message_string = 'dissipation_1d = "' // TRIM( dissipation_1d ) //      &
3569                        '" is unknown'
3570       CALL message( 'check_parameters', 'PA0138', 1, 2, 0, 6, 0 )
3571    ENDIF
3572
3573!
3574!-- Set time for the next user defined restart (time_restart is the
3575!-- internal parameter for steering restart events)
3576    IF ( restart_time /= 9999999.9_wp )  THEN
3577       IF ( restart_time > time_since_reference_point )  THEN
3578          time_restart = restart_time
3579       ENDIF
3580    ELSE
3581!
3582!--    In case of a restart run, set internal parameter to default (no restart)
3583!--    if the NAMELIST-parameter restart_time is omitted
3584       time_restart = 9999999.9_wp
3585    ENDIF
3586
3587!
3588!-- Set default value of the time needed to terminate a model run
3589    IF ( termination_time_needed == -1.0_wp )  THEN
3590       IF ( host(1:3) == 'ibm' )  THEN
3591          termination_time_needed = 300.0_wp
3592       ELSE
3593          termination_time_needed = 35.0_wp
3594       ENDIF
3595    ENDIF
3596
3597!
3598!-- Check the time needed to terminate a model run
3599    IF ( host(1:3) == 't3e' )  THEN
3600!
3601!--    Time needed must be at least 30 seconds on all CRAY machines, because
3602!--    MPP_TREMAIN gives the remaining CPU time only in steps of 30 seconds
3603       IF ( termination_time_needed <= 30.0_wp )  THEN
3604          WRITE( message_string, * )  'termination_time_needed = ',            &
3605                 termination_time_needed, ' must be > 30.0 on host "',         &
3606                 TRIM( host ), '"'
3607          CALL message( 'check_parameters', 'PA0139', 1, 2, 0, 6, 0 )
3608       ENDIF
3609    ELSEIF ( host(1:3) == 'ibm' )  THEN
3610!
3611!--    On IBM-regatta machines the time should be at least 300 seconds,
3612!--    because the job time consumed before executing palm (for compiling,
3613!--    copying of files, etc.) has to be regarded
3614       IF ( termination_time_needed < 300.0_wp )  THEN
3615          WRITE( message_string, * )  'termination_time_needed = ',            &
3616                 termination_time_needed, ' should be >= 300.0 on host "',     &
3617                 TRIM( host ), '"'
3618          CALL message( 'check_parameters', 'PA0140', 1, 2, 0, 6, 0 )
3619       ENDIF
3620    ENDIF
3621
3622!
3623!-- Check pressure gradient conditions
3624    IF ( dp_external .AND. conserve_volume_flow )  THEN
3625       WRITE( message_string, * )  'Both dp_external and conserve_volume_flo', &
3626            'w are .TRUE. but one of them must be .FALSE.'
3627       CALL message( 'check_parameters', 'PA0150', 1, 2, 0, 6, 0 )
3628    ENDIF
3629    IF ( dp_external )  THEN
3630       IF ( dp_level_b < zu(nzb) .OR. dp_level_b > zu(nzt) )  THEN
3631          WRITE( message_string, * )  'dp_level_b = ', dp_level_b, ' is out ', &
3632               ' of range'
3633          CALL message( 'check_parameters', 'PA0151', 1, 2, 0, 6, 0 )
3634       ENDIF
3635       IF ( .NOT. ANY( dpdxy /= 0.0_wp ) )  THEN
3636          WRITE( message_string, * )  'dp_external is .TRUE. but dpdxy is ze', &
3637               'ro, i.e. the external pressure gradient & will not be applied'
3638          CALL message( 'check_parameters', 'PA0152', 0, 1, 0, 6, 0 )
3639       ENDIF
3640    ENDIF
3641    IF ( ANY( dpdxy /= 0.0_wp ) .AND. .NOT. dp_external )  THEN
3642       WRITE( message_string, * )  'dpdxy is nonzero but dp_external is ',     &
3643            '.FALSE., i.e. the external pressure gradient & will not be applied'
3644       CALL message( 'check_parameters', 'PA0153', 0, 1, 0, 6, 0 )
3645    ENDIF
3646    IF ( conserve_volume_flow )  THEN
3647       IF ( TRIM( conserve_volume_flow_mode ) == 'default' )  THEN
3648
3649          conserve_volume_flow_mode = 'initial_profiles'
3650
3651       ELSEIF ( TRIM( conserve_volume_flow_mode ) /= 'initial_profiles' .AND.  &
3652            TRIM( conserve_volume_flow_mode ) /= 'inflow_profile' .AND.        &
3653            TRIM( conserve_volume_flow_mode ) /= 'bulk_velocity' )  THEN
3654          WRITE( message_string, * )  'unknown conserve_volume_flow_mode: ',   &
3655               conserve_volume_flow_mode
3656          CALL message( 'check_parameters', 'PA0154', 1, 2, 0, 6, 0 )
3657       ENDIF
3658       IF ( (bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic')  .AND.                &
3659          TRIM( conserve_volume_flow_mode ) == 'bulk_velocity' )  THEN
3660          WRITE( message_string, * )  'non-cyclic boundary conditions ',       &
3661               'require  conserve_volume_flow_mode = ''initial_profiles'''
3662          CALL message( 'check_parameters', 'PA0155', 1, 2, 0, 6, 0 )
3663       ENDIF
3664       IF ( bc_lr == 'cyclic'  .AND.  bc_ns == 'cyclic'  .AND.                 &
3665            TRIM( conserve_volume_flow_mode ) == 'inflow_profile' )  THEN
3666          WRITE( message_string, * )  'cyclic boundary conditions ',           &
3667               'require conserve_volume_flow_mode = ''initial_profiles''',     &
3668               ' or ''bulk_velocity'''
3669          CALL message( 'check_parameters', 'PA0156', 1, 2, 0, 6, 0 )
3670       ENDIF
3671    ENDIF
3672    IF ( ( u_bulk /= 0.0_wp .OR. v_bulk /= 0.0_wp ) .AND.                      &
3673         ( .NOT. conserve_volume_flow .OR.                                     &
3674         TRIM( conserve_volume_flow_mode ) /= 'bulk_velocity' ) )  THEN
3675       WRITE( message_string, * )  'nonzero bulk velocity requires ',          &
3676            'conserve_volume_flow = .T. and ',                                 &
3677            'conserve_volume_flow_mode = ''bulk_velocity'''
3678       CALL message( 'check_parameters', 'PA0157', 1, 2, 0, 6, 0 )
3679    ENDIF
3680
3681!
3682!-- Check particle attributes
3683    IF ( particle_color /= 'none' )  THEN
3684       IF ( particle_color /= 'absuv'  .AND.  particle_color /= 'pt*'  .AND.   &
3685            particle_color /= 'z' )  THEN
3686          message_string = 'illegal value for parameter particle_color: ' //   &
3687                           TRIM( particle_color)
3688          CALL message( 'check_parameters', 'PA0313', 1, 2, 0, 6, 0 )
3689       ELSE
3690          IF ( color_interval(2) <= color_interval(1) )  THEN
3691             message_string = 'color_interval(2) <= color_interval(1)'
3692             CALL message( 'check_parameters', 'PA0315', 1, 2, 0, 6, 0 )
3693          ENDIF
3694       ENDIF
3695    ENDIF
3696
3697    IF ( particle_dvrpsize /= 'none' )  THEN
3698       IF ( particle_dvrpsize /= 'absw' )  THEN
3699          message_string = 'illegal value for parameter particle_dvrpsize:' // &
3700                           ' ' // TRIM( particle_color)
3701          CALL message( 'check_parameters', 'PA0314', 1, 2, 0, 6, 0 )
3702       ELSE
3703          IF ( dvrpsize_interval(2) <= dvrpsize_interval(1) )  THEN
3704             message_string = 'dvrpsize_interval(2) <= dvrpsize_interval(1)'
3705             CALL message( 'check_parameters', 'PA0316', 1, 2, 0, 6, 0 )
3706          ENDIF
3707       ENDIF
3708    ENDIF
3709
3710!
3711!-- Check nudging and large scale forcing from external file
3712    IF ( nudging .AND. ( .NOT. large_scale_forcing ) )  THEN
3713       message_string = 'Nudging requires large_scale_forcing = .T.. &'//      &
3714                        'Surface fluxes and geostrophic wind should be &'//    &
3715                        'prescribed in file LSF_DATA'
3716       CALL message( 'check_parameters', 'PA0374', 1, 2, 0, 6, 0 )
3717    ENDIF
3718
3719    IF ( large_scale_forcing .AND. ( bc_lr /= 'cyclic'  .OR.                   &
3720                                    bc_ns /= 'cyclic' ) )  THEN
3721       message_string = 'Non-cyclic lateral boundaries do not allow for &' //  &
3722                        'the usage of large scale forcing from external file.'
3723       CALL message( 'check_parameters', 'PA0375', 1, 2, 0, 6, 0 )
3724     ENDIF
3725
3726    IF ( large_scale_forcing .AND. ( .NOT. humidity ) )  THEN
3727       message_string = 'The usage of large scale forcing from external &'//   & 
3728                        'file LSF_DATA requires humidity = .T..'
3729       CALL message( 'check_parameters', 'PA0376', 1, 2, 0, 6, 0 )
3730     ENDIF
3731
3732    IF ( large_scale_forcing .AND. topography /= 'flat' )  THEN
3733       message_string = 'The usage of large scale forcing from external &'//   & 
3734                        'file LSF_DATA is not implemented for non-flat topography'
3735       CALL message( 'check_parameters', 'PA0377', 1, 2, 0, 6, 0 )
3736    ENDIF
3737
3738    IF ( large_scale_forcing .AND.  ocean  )  THEN
3739       message_string = 'The usage of large scale forcing from external &'//   & 
3740                        'file LSF_DATA is not implemented for ocean runs'
3741       CALL message( 'check_parameters', 'PA0378', 1, 2, 0, 6, 0 )
3742    ENDIF
3743
3744    CALL location_message( 'finished', .TRUE. )
3745
3746!
3747!-- Prevent empty time records in volume, cross-section and masked data in case of
3748!-- non-parallel netcdf-output in restart runs
3749    IF ( netcdf_data_format < 5 )  THEN
3750       IF ( TRIM( initializing_actions ) == 'read_restart_data' )  THEN
3751          do3d_time_count    = 0
3752          do2d_xy_time_count = 0
3753          do2d_xz_time_count = 0
3754          do2d_yz_time_count = 0
3755          domask_time_count  = 0
3756       ENDIF
3757    ENDIF
3758
3759!
3760!-- Check &userpar parameters
3761    CALL user_check_parameters
3762
3763
3764 END SUBROUTINE check_parameters
Note: See TracBrowser for help on using the repository browser.