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

Last change on this file since 1554 was 1554, checked in by maronga, 9 years ago

last commit documented

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