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

Last change on this file since 1557 was 1557, checked in by suehring, 9 years ago

Enable monotone advection for scalars in combination with fifth-order scheme using monotonic limiter.

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