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

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

bugfix: compilation of radiation_model.f90 without netcdf

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