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

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

last commit documented

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