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

Last change on this file since 1560 was 1558, checked in by suehring, 10 years ago

last commit documented

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