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

Last change on this file since 1763 was 1763, checked in by hellstea, 8 years ago

last commit documented

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