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

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

Added support for RRTMG radiation code

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