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

Last change on this file since 1748 was 1746, checked in by gronemeier, 8 years ago

last commit documented

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